zyndagj / BSMAPz

Updated and optimized fork of BSMAP
Other
18 stars 6 forks source link

Questions about the bam output #4

Closed dingailuma closed 6 years ago

dingailuma commented 6 years ago

Hi, Recently I've tried using bsmap to align bisulfite sequencing data, I found something I can't understand when looking at the output bam file:

ST-E00206:158:HVLMYCCXX:7:121   161 # chrM  156 255 106M    # chr18 29124951    0   ATTTATCACACCTACATTCAATATTACAAACAAACATACCTACTAAAATATATTAATTAATTAATACTTATAAAACATAATAATAACAATTAAATATCTACACAAC  AA<FFAAAFKKFAFKKFKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK  NM:i:0  ZS:Z:--
ST-E00206:158:HVLMYCCXX:7:1211  129 # chrM  156 255 144M    # chr1  161734142   0   ATTTATCACAGCTACAGTCAATGTTACGAACAGACATACTTACTAAAATATATTAATTAATTAATACTTATGAAACATAATAATAACAATTAAATATCTACACAACCACTTTCCACACAAACATCATAACAAAAAATTTCCACC    A<FFFFF,FAFA,AFF,F<F<K<FKK,,F,,F,,,A<,,,A<,7,7AA7FAA,FA,FK,AF7,7,,7,F,F,,F,AF<FKKKAA<FKKFAFFFAF,FAAAFKKAFFKKKKFKFFAKFFKF7F7AAFFKKKKKFKKKKKF7AAFA    NM:i:7  ZS:Z:--
ST-E00206:158:HVLMYCCXX:7:1211  129 # chrM  156 255 144M    # chr18 1125144 0   AGTTATCACACCTACATTCAATATTACAAACAAACATACTTACTAAAATATATTAATTAATTAATACTTATAAAACATAATAATAACAATTAAATATCTACACAACCACTTTCCACACAAACATCATAACAAAACATTTCCACC    A<AFAA,F,,FAAFKKKF<<,FFF<A,<<FAAFFFFKFFKK<,FFFKKKKKFKKKKKKKKKKKKKFKFK7AKA7AFKFFFKFKKKKKKKKKFAKKKKFAKAFKKKFFKFFFKKAKFAFFKFKKKKKKFKKKKFK,FKKKFFAAA    NM:i:3  ZS:Z:--

Above is from the bam file cut, could you please explain why in the same line of bam, different chromosomes appear ? In my understanding, read1 and read2 should be on the same chromosome, is that right ? Thank you so much !

zyndagj commented 6 years ago

Reads are usually aligned independently (unpaired) and then filtered downstream. You can discard pairs like this with either samtools

samtools view -bf 0x2 in.bam > properly_paired.bam

or even using methratio.py from this repository

methratio.py -p in.bam > out.txt

You should also consider filtering out non-unique alignments. You can read more about sam/bam flags at http://www.samformat.info/sam-format-flag or the main samtools documentation page.

dingailuma commented 6 years ago

Thank you so much !