The -F 2820 bit might not be necessary depending on the alignment, these could be primary mapped reads. I know when I process my exomes I map against the whole genome because of the way the baited capture fragments work. I expect to get primary, proper alignments that would meet these criteria.
extract reads where one pair of read is at the annotated region and other not
Hello!
I have a bam alignment file from lifescope. I would like to extract reads that fulfil following criteria:
- one pair is aligned to the annotated region (on exom) and
- the other pair aligned to the intron region
I tried the following command but was unable to grep the condition.
samtools view -F 14 input_file.bam | grep -v " = " > ouput_file.sam
I have searched for the answer and played with different options for samtools view but failed to find that satisfies the criteria.
Thanks
• 2,677 views
•
link
updated
by
Ram
•
written
by
nikkihathi •
•
written
by
nikkihathi •
1 answer
get the read names mapping the exons, and sort
samtools view -F 2820 -L exon.bed in.bam | cut -f 1 | LC_ALL=C sort > name1.txt
get the read names mapping the introns, and sort
samtools view -F 2820 -L intron.bed in.bam | cut -f 1 | LC_ALL=C sort > name2.txt
2820 : read unmapped, not primary alignment , read fails platform/vendor quality checks supplementary alignment
get the common read names:
comm -12 name1.txt name2.txt
go back to the bam, and get the reads with those names.
EDIT: you might get some reads overlapping BOTH exon and intron. you can play with bedtools to remove those reads...
• 1 views
•
link
written
by
Pierre Lindenbaum
Log in to answer this question.
More posts like this
-
Filtering a 10X generated .bam file based on a list of barcodes
written by jhnicolas •
-
Samtools merge bam issue with header tags
written by sam •
-
How To Extract Reads-Pairs Aligned Concordantly Exactly 1 Time from Bowtie2 (version 2.3.4.1 or hig…
written by Chilly •
-
samtools view region and include duplicate reads
written by mgene •
-
low alignment rate with bowtie2 of ChIP-seq data 150bp pair end
written by Shuang He •
-
Extract only RNA-seq reads in a BAM file that extend into an intron
written by foxdie •
-
Using samtools to view an indexed file in a private S3 bucket
written by Brynjar •
-
differential RNA seq analysis
written by chrisclarkson100
-
Samtools view command, Exec format error while trying to view a human chromosome region
written by ricfoz •
-
Coverage depth and breadth after recaliberation
written by nikkihathi •
There's no possible way to do this with a combination of samtools and grep. You're going to write a script using pysam or something similar to do this.