Hot answers tagged

18

Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via #!/bin/sh or e.g. #!/bin/bash if it matters), so won't be affected by your choice of user shell. If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell. See Csh Programming Considered Harmful and ...


11

TL;DR SH adheres to an official industry standard, but it is not suitable for scientific computing. Bash is considered an informal standard (e.g., by Google). Bash 3 is preferable in most of situations in the world of bioinformatics. Long answer As already described in other answers, SH (/bin/sh, plain Bourne shell, the original UNIX shell) should fully ...


9

Here's a simple awk approach: awk '{if(/^>/){split($1,a,"[|.]")}print >> a[2]".fa"}' Protein_FASTA.txt Or, more concisely, just: awk '/^>/{split($1,a,"[|.]")}{print >> a[2]".fa"}' Protein_FASTA.txt When run on the file linked to in your question, that results in the following files: $ ls AP014650.fa CP003848.fa CP007713.fa ...


9

Use perl-rename. This is usually called rename on Debian-based systems like Ubuntu or Mint, and perl-rename or prename on others. Assuming you have it as rename, you can simply do: rename -n 's/.*\.vcf/"A" . ++$c . ".vcf"/e' *snp.pass.vcf The -n causes rename to only print what it would do, without doing anything. So once you make sure it does what you ...


8

All you need is cat. You won't find any better tool for a simple job like this. Just run: cat NA24694_GCCAAT_L001_R1*fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz cat NA24694_GCCAAT_L001_R2*fastq.gz > EA00694_GCCAAT_L001_R2.fastq.gz cat NA24694_GCCAAT_L001_R3*fastq.gz > EA00694_GCCAAT_L002_R1.fastq.gz cat NA24694_GCCAAT_L001_R4*fastq.gz > ...


7

Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve. Since I don’t know how how to do this in sed, and it’s not trivial in awk either, here’s a Perl ...


7

The Open Group Base Specifications Issue 7 IEEE Std 1003.1™-2008, 2016 Edition, or "The POSIX Standard" for short, is the standard that defines the interfaces and utilities provided by a Unix system. Among these is the command line shell language and tools (see "Shell & Utilities" in the main index on the page linked above). As far as I know, there is ...


6

Just today I wrote a script to do exactly this using Biopython. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script filter_fasta_by_list_of_headers.py is : #!/usr/bin/env python3 from Bio import SeqIO import sys ffile = SeqIO.parse(sys.argv[1], "fasta") header_set = set(line.strip() for ...


6

Suppose you keep sequence names in ids.txt and sequences in seq.fa: awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa


6

If you want to learn how to do things with command-line tools, you can linearize the FASTA with awk, pipe to grep to filter for items of interest named in patterns.txt, then pipe to tr to delinearize: $ awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff ...


6

You have an extra done, instead you want: for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq); do echo $i bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}_R1.fq.gz | samtools view -bS - >${i}.bam done As an aside, I'm not a big fan of piping ls to cut and then uniq. Why not instead: for f in *.fq.gz; do sampleName=${f%%_R1.fq....


5

Your problem is that you have written a script in bash but are then running it using sh. Bash, the Bourne-again shell, is not the same as the Bourne shell (sh). What is more, on Ubuntu, sh is actually not even the Bourne shell but another minimal shell called dash. The /bin/sh on Ubuntu is a symlink to dash: $ ls -l /bin/sh lrwxrwxrwx 1 root root 4 Mar 1 ...


5

The generic command sh is quite literally an industry standard, a POSIX standard, to be precise (IEEE 1003.2 and 1003.2a, available for purchase for hundreds of dollars at various websites). In theory, any script that starts with #!/bin/sh should conform to this standard. In practise, most Linux systems have a shell that is close to this standard, but has a ...


5

I would not say bash as a "standard", but it is indeed likely to be the most widely used unix shell and available by default on most modern unix/linux distros. There are a few other more convenient shells like zsh that are broadly compatible with /bin/sh, but they are not as widely available. There is also C-shell and in particular its open-source ...


5

Perl should be fairly fast with this when using a hash set to store the list of lines. A structure like this also works for subsetting based on a field value, where the comparison would be with the field rather than "$.": #!/usr/bin/perl use strict; use warnings; my $lines_file = $ARGV[0]; my %include_lines = (); open my $lines_fh, '<', $lines_file or ...


5

just to complete your question, you can do it in R in one line. After setting the directory as your working directory, just type the following: file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf")) Maybe you can try file.rename(list.files(pattern = ".vcf"), ...


4

The following sed command will do the trick. sed -e '/^[^>]/ s/[a-z]/\-/g' in.fna > out.fna Lines starting with > will be escaped and the lowercase characters will be replaced with - .


4

Some related questions appear in other sites, with potentially interesting solutions, which I report here: To sample approximately 1% of the non-empty lines: awk 'BEGIN {srand()} !/^$/ { if (rand() <= .01) print $0}' input_file (from https://stackoverflow.com/a/692321/1878788) To select 1000 random lines: shuf -n 1000 input_file (from https://...


4

echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP 'TTCG[CGT][ATGC]{1,15}TTCG[CGT]' [CGT] looks for one occurrence of one of the three. [ATGC]{1,15} looks for "up to 15 letter long" combinations of the four bases. Moreover, -P makes your regex call Perl-like and is required to make this regex work. Different regex flavors ...


4

$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f. EDIT AS THE OP HAS EDITED THE QUESTION: First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the ...


3

A very simple approach, but one which will work with the names you show, is to remove any occurrence of 00N (where N is any number) from right before the .gz. To do so, edit your script to: for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \ -1 ${i}...


3

We sequence DNA on daily basis. Like terdon said, you can use cat to group all the FastQ of all your sample whether they are bgzipped or not like so: root@slurm1: ~ # ls am224_S6_L001_R1_001.fastq.gz am224_S6_L001_R2_002.fastq.gz am224_S6_L002_R1_003.fastq.gz am224_S6_L002_R2_004.fastq.gz am232_S10_L001_R1_001.fastq.gz am232_S10_L001_R2_002.fastq.gz ...


3

If I understand correctly, your question boils down to: parsing a bunch of numbers in a row from a .csv file gives a different result from parsing a bunch of numbers as rows of a text file. I that right? I suggest that you look at the source of your text file. For example, was the text file created on a Windows host (running bash, or in the command ...


3

You could try it with sed. Replace all lower-case chars with - in lines not starting with >: sed -e '/^>/!s/[a-z]/-/g' in.fa


3

Using the FastaToTbl and TblToFasta scripts I have posted before, you can do: FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta > file.2.fa


3

If your filenames do not contain blank spaces, you can do it with a for-loop in bash: Make a test run with echo so that the command is only printed but not yet executed: for i in *.vcf; do echo mv $i A${k}.vcf ; let k++ ; done And if the result looks good remove the echo statement to execute the mv commands: for i in *.vcf; do mv $i A${k}.vcf ; let k++ ; ...


3

EDIT: (Based on the comments by @terdon, with minor changes) Try this Perl one-liner. It is recommended for general use, even if your file names contain "unexpected" characters such as newlines. perl -e 'for ( 0..$#ARGV ) { rename $ARGV[$_] => "A@{[ $_ + 1 ]}.vcf"; } ' *pass.vcf Here, -e tells perl to look for code on the command line instead of a ...


3

You don't need a loop. You can do the whole thing with a simple awk one-liner: awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta I ran this on your files and got: $ for file in gene*; do echo "=== File: $file ==="; cat $file; done === File: gene1.fasta === >locus1 ATGCGTAGAG >ysr ATGTAGCGA >siv TAGTAGTAT === ...


Only top voted, non community-wiki answers of a minimum length are eligible