bowtie2 output only valid alignments
1
0
Entering edit mode
8.0 years ago
ddzhangzz ▴ 90

I was trying to align paired end RNASeq data to a locally-built reference using bowtie2:

$bowtie2 --trim3 0 --threads 6 --no-discordant --no-mixed -k 1 \
    -x $idx -1 $f1 -2 $f2 -S classI.sam

And the sam file output looks like this (ignored the header):

ERR356371.7445  141 *   0   0   *   *   0   0   ATTGAGAAGCCCCCCTTCGAGCTGCCAGACTTCATCAAACGCACAGGCA   2@@@DDDDA?DFDHGBHHICEGGDFDGCCGIIIJJBFHD<DFHIJAAGH   YT:Z:UP
ERR356371.7446  77  *   0   0   *   *   0   0   CTACCATCACCCTCTACATCACCGCCCCGACCTTAGCTCTCACCATCGC   CCCFFFFFHHHHHJJIJJJJIJJIJJJJJIJJJJJJJJIJGGIJEHJJI   YT:Z:UP
ERR356371.7446  141 *   0   0   *   *   0   0   ATCGCAGTGCGCCGATCAGGGCGTAGTTTGAGTTTGATGCTCACCCTGA   =CCCFFFDFHHHHHGIJJJJJJJHIGFIIJGIIIJJIJJJJJJJJJJJH   YT:Z:UP
.
.
ERR356371.7400  99  HLA:HLA05978    527 255 49M =   636 158 CGGCGGAGCAGCGGAGAGTCTACCTGGAGGGCCGGTGCGTGGACGGGCT   CCCFFFFFGHHGHJJIIJHFHJIIJJJJJJJJJHF=BDDBDDDDDDDDD   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YS:i:-8 YT:Z:CP
ERR356371.7400  147 HLA:HLA05978    636 255 49M =   527 -158    TATGACCCACCACCCCATCTCTGACCATGAGGCCACCCTGAGGTGCAGT   JJIHGJHF@HD?JJIJIHHGJIIGJJJIJJJJIGEHHHHHFFFFFCCC2   AS:i:-8 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:46T1G0 YS:i:0  YT:Z:CP
.
.
.
ERR356371.7447  77  *   0   0   *   *   0   0   AATTCATGTCTGCACCAACTTGGACAACTGGGGGAGGCTGGATGCATCA   CCCFFFFFHHHHHJJJJJJJJJJIJJJJJJJJJJGHJIJIJGJJJIJJI   YT:Z:UP
ERR356371.7447  141 *   0   0   *   *   0   0   ATTGAAGTCAAGGTGTTGGCTGCCACATGCCATGTCCTCTTAACAAGAT   =CCCFFFFFHHHHFGIJJJJJJJJIJJJJJJIJJIJJJJJJJJJJJJJJ   YT:Z:UP
ERR356371.7448  77  *   0   0   *   *   0   0   AATTGCCAATATTTGCCCCTCCCCTTTAATAAAACTTTTAAGAAGTCAC   <?@DDFFDHFFFHJGEFBFFGAEDGHGIJIFIEF@GHGHBHEBFHGFIH   YT:Z:UP
ERR356371.7448  141 *   0   0   *   *   0   0   AGGCAAATATTGGCAATTAGTTGGCCGTAAGTACAGCACAGTGCAGCTT   A??@D?DDDDFHDH@GGGIGAACDBHGCCEFIIIIG@DDDBHGGGGGHI   YT:Z:UP

My question is whether there is an option in bowtie2 only output the valid alignments such as:

ERR356371.7400  99  HLA:HLA05978    527 255 49M =   636 158 CGGCGGAGCAGCGGAGAGTCTACCTGGAGGGCCGGTGCGTGGACGGGCT   CCCFFFFFGHHGHJJIIJHFHJIIJJJJJJJJJHF=BDDBDDDDDDDDD   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YS:i:-8 YT:Z:CP
ERR356371.7400  147 HLA:HLA05978    636 255 49M =   527 -158    TATGACCCACCACCCCATCTCTGACCATGAGGCCACCCTGAGGTGCAGT   JJIHGJHF@HD?JJIJIHHGJIIGJJJIJJJJIGEHHHHHFFFFFCCC2   AS:i:-8 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:46T1G0 YS:i:0  YT:Z:CP

Namely the middle part of the full results. Does somebody know how to output this?

RNA-Seq • 2.6k views
ADD COMMENT
3
Entering edit mode
8.0 years ago
mastal511 ★ 2.1k

There seems to be an option in the manual, under the list of parameters, and 'SAM options',

--no-unal    Suppress SAM records for reads that failed to align.

That may do what you want.

ADD COMMENT
1
Entering edit mode

it does make more sense not to print out unaligned reads if you already know that you're not interested in them, plus it would be faster. in case you aren't completely sure about it, I would suggest to convert the sam file into a bam file (samtools sort followed by samtools view -Sb) and then extract from that bam file only the aligned reads (samtools view -F4). you will end up having to deal with temporal files though.

ADD REPLY

Login before adding your answer.

Traffic: 1924 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6