How to extra read pair orientation from paired end reads

Question

[Wang Cong]

Hi, I have bam files for paired-end sequencing from the Hi-C experiment. I want to know the fraction of read pair orientation as a function of mapping distance. May I ask if there is software for this?

If not, I first need to extract different read pair orientations (forward-reverse, forward-forward, reverse-forward, reverse-reverse). May I ask if is a way to do this? Thanks.

Answer 1

[Pierre Lindenbaum]

I want to know the fraction of read pair orientation as a function of mapping distance.

that's not clear to me.

If not, I first need to extract different read pair orientations (forward-reverse, forward-toward, reverse-toward, reverse-reverse).

use samtools view with a FILTER EXPRESSION.

 samtools view --expr 'flag.paired && !flag.unmap && !flag.munmap && flag.reverse && !flag.mreverse' in.bam

[Wang Cong]

Thanks! So this command extracts reverse-forward pairs?

[Pierre Lindenbaum]

Thanks! So this command extracts reverse-forward pairs?

reverse R1, forward R2

[Wang Cong]

Thank you! I used the command. However, I found the output reads are not paired. The inputs are all paired. All reads are supposed to be paired, I think?

Output:

58201619 + 0 in total (QC-passed reads + QC-failed reads)
58201619 + 0 primary
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
58201619 + 0 mapped (100.00% : N/A)
58201619 + 0 primary mapped (100.00% : N/A)
58201619 + 0 paired in sequencing
29103737 + 0 read1
29097882 + 0 read2
0 + 0 properly paired (0.00% : N/A)
58201619 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
10742519 + 0 with mate mapped to a different chr
10742519 + 0 with mate mapped to a different chr (mapQ>=5)

Input:

116394722 + 0 in total (QC-passed reads + QC-failed reads)
116394722 + 0 primary
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
116394722 + 0 mapped (100.00% : N/A)
116394722 + 0 primary mapped (100.00% : N/A)
116394722 + 0 paired in sequencing
58197361 + 0 read1
58197361 + 0 read2
0 + 0 properly paired (0.00% : N/A)
116394722 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
21484650 + 0 with mate mapped to a different chr
21484650 + 0 with mate mapped to a different chr (mapQ>=5)

[Pierre Lindenbaum]

I found the output reads are not paired.

uh ??

58201619 + 0 in total (QC-passed reads + QC-failed reads)
58201619 + 0 paired in sequencing

[Wang Cong]

Sorry, I think the read number should be even since all the reads are in pairs but now it's an odd number. I think read 1 and read 2 should be discarded together according to the command?