Calculate the number of paired-end reads their middle size is greater than # KB

Question

[lisadavic66]

Dear All,

I have a BAM file of paired-end sequencing reads. I want to calculate how many paired-end reads that mapped to the same chromosomes their middle size (the number of base pairs) between the two ends of paired-end reads is greater than 100 kb, and how many paired-end reads that mapped to the different chromosomes their middle size (the number of base pairs) between the two ends of paired-end reads is greater than 100 kb. Can anybody help me get this done? Thank you very much in advance.

Answer 1

[James Ashmore]

You can use a program called BEDTools to do this analysis. First convert your BAM to a BAMPE file, then use AWK to get the number of bases between the two pairs with a cutoff. For example:

$ bedtools bamtobed -i reads.bam -bedpe > reads.bedpe
chr7   118965072   118965122   chr7   118970079   118970129 TUPAC_0001:3:1:0:1452#0 37     +     -
$ awk -v CUTOFF=500 -v OFS="\t" '$5 - $2 >= CUTOFF {print $1, $5 - $2}' reads.bedpe > distances.txt
chr7 5007

[lisadavic66]

Thank you so much for your help. I will try your script.

How can I calculate the numbers for all of chromosomes efficiently? and how about that mapped to different chromosomes? Thank you again.

[James Ashmore]

You can cut the chromosome field out from the distances.txt file and then count how often each chromosome identifier appears. For example:

cut -f1 distances.txt | sort | uniq -c

[Antonio R. Franco]

You can easily discover the number of reads whose mate mapped to a different chromosome with

samtools file.bam flagstat