added 1 character in body
Source Link
terdon
  • 7.8k
  • 3
  • 14
  • 40

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (`rm *tmprm *tmp).

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (`rm *tmp).

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (rm *tmp).

added 633 characters in body
Source Link
terdon
  • 7.8k
  • 3
  • 14
  • 40

As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.

Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (>chrN) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.

But I don't really know why you would want to do it. You're the one asking for it, after all! :)


As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.

Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (>chrN) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.

But I don't really know why you would want to do it. You're the one asking for it, after all! :)

added 2073 characters in body
Source Link
terdon
  • 7.8k
  • 3
  • 14
  • 40

Using the TblToFasta and FastaToTbl scriptscripts I have posted before, combined with classic *nix utils, you can do:

###Explanation

The FastaToTbl will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to join, which will join its input on lines where the 1st field is the same (by default). The -a1 -a2 means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.

Since the join will add a space between the two joined lines, we deal with that by replacing the second space on each line with a = (sed 's/ /=/2') and then pass the whole thing to TblToFasta to revert to fasta format. Finally, we use tr to replace the = with a newline so you can see where each sequence came from as you show in your example.

###Explanation Here, the tbl format sequences (see above) are passed to an awk script which does the work:

  • seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2 : if the 1st field ($1, the sequence identifier) is already in the array seq, then set the value of seq[$1] to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field ($2, the sequence). If it hasn't been seen before, then set the value of seq[$1] to the 2nd field, the sequence.
  • END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]} : after processing everything, go through each element in the array seq and print out the data in fasta format.

###Explanation This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

  • if(!seen[$1]++){print ">"$1 > $1".tmp"} : if this is the first time we've seen this ID, print it as a fasta header into the file $1".tmp" (e.g. chr1.tmp).
  • print $2 >> $1".tmp" append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (`rm *tmp).

Using the TblToFasta and FastaToTbl script I have posted before, combined with classic *nix utils, you can do:

Using the TblToFasta and FastaToTbl scripts I have posted before, combined with classic *nix utils, you can do:

###Explanation

The FastaToTbl will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to join, which will join its input on lines where the 1st field is the same (by default). The -a1 -a2 means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.

Since the join will add a space between the two joined lines, we deal with that by replacing the second space on each line with a = (sed 's/ /=/2') and then pass the whole thing to TblToFasta to revert to fasta format. Finally, we use tr to replace the = with a newline so you can see where each sequence came from as you show in your example.

###Explanation Here, the tbl format sequences (see above) are passed to an awk script which does the work:

  • seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2 : if the 1st field ($1, the sequence identifier) is already in the array seq, then set the value of seq[$1] to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field ($2, the sequence). If it hasn't been seen before, then set the value of seq[$1] to the 2nd field, the sequence.
  • END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]} : after processing everything, go through each element in the array seq and print out the data in fasta format.

###Explanation This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

  • if(!seen[$1]++){print ">"$1 > $1".tmp"} : if this is the first time we've seen this ID, print it as a fasta header into the file $1".tmp" (e.g. chr1.tmp).
  • print $2 >> $1".tmp" append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (`rm *tmp).

Source Link
terdon
  • 7.8k
  • 3
  • 14
  • 40
Loading