Hi there,
I've been stumbling through some Illumina HiSeq data for a while, and while I've managed to answer some of the questions we wanted to ask, I'm having trouble with some others. One such question is whether it is possible to identify and print out only those paired-end reads that map to different chromosomes or have a larger insert size than expected. I was hoping to use this for identifying new retrotransposition events and rare fusions.
I have aligned the reads using bwa sampe, and aligned the forward and reverse separately with bwa samse.
My questions are as follows:
1. Can I use my paired-end aligned sam/bam file to pull out the reads where each mate maps to a different chromosome? I'm not sure if bwa biases against these and forces poorer alignments that are spatially close (I know that ambiguously aligned reads have the XA:Z: column that identifies all other possible mapping locations, but not sure if I can use this). If the sam file I have can't be used to answer this question, could anyone suggest a good program (fusionhunter?) that can do this?
2. If I can use the sam file I currently have (and not have to re-align with a specialized program), could anyone suggest basic code to pull out these sequences? I'm thinking of using either awk, Perl or R, but I'm still learning how to write in these languages.
I would want to match lines that have the same first column of the sam file (ignoring the #0 and #0/3 at the end which denotes forward or reverse read), then only pull out the lines where the thrid column (chromosome) is different (or the forth column (location in bp) is greater than 50kb).
Any help or advice would be greatly appreciated.
Thanks fo your time.
I've been stumbling through some Illumina HiSeq data for a while, and while I've managed to answer some of the questions we wanted to ask, I'm having trouble with some others. One such question is whether it is possible to identify and print out only those paired-end reads that map to different chromosomes or have a larger insert size than expected. I was hoping to use this for identifying new retrotransposition events and rare fusions.
I have aligned the reads using bwa sampe, and aligned the forward and reverse separately with bwa samse.
My questions are as follows:
1. Can I use my paired-end aligned sam/bam file to pull out the reads where each mate maps to a different chromosome? I'm not sure if bwa biases against these and forces poorer alignments that are spatially close (I know that ambiguously aligned reads have the XA:Z: column that identifies all other possible mapping locations, but not sure if I can use this). If the sam file I have can't be used to answer this question, could anyone suggest a good program (fusionhunter?) that can do this?
2. If I can use the sam file I currently have (and not have to re-align with a specialized program), could anyone suggest basic code to pull out these sequences? I'm thinking of using either awk, Perl or R, but I'm still learning how to write in these languages.
I would want to match lines that have the same first column of the sam file (ignoring the #0 and #0/3 at the end which denotes forward or reverse read), then only pull out the lines where the thrid column (chromosome) is different (or the forth column (location in bp) is greater than 50kb).
Any help or advice would be greatly appreciated.
Thanks fo your time.
Comment