Hello everybody,
it is possible to get only the reads if the start coordinates is in a region.
When i run this: samtools view cigarN.bam X:68206675-68208264
I get this:
FCC0WRYACXX:1:1103:14931:24492 99 X 68204511 255 59M2530N30M = 68207133 2703 TGTTCTTATTATTATGGTATCATTTATTTGTAAAATTTATAATTGTAAAATACAAACTAGTAAGAACCCTATCAGTGTTGTTTTTGTTT HHJJJJJJJJJJJIJJJFHIJJJJJIJJJJJIJJJJJJJJJJJIJJJJJJJJJJJJJIJJIIIJGGIIIJJJJJJEHFHHHFFFFDEEE PG:Z:MarkDuplicates NH:i:1 HI:i:1 nM:i:1 AS:i:165
FCC0WRYACXX:1:1304:17536:3030 99 X 68206627 255 47M3479N42M = 68210195 3656 TTTCTCATACTGGCAATAGTATTTCACCCAGGATATTCCTAAAGCCCATTTCTCTTGTGTATAGAGATAGCCTTCAATAGTTGGCTGTC HHJJJJIJJJJJJJJIIIJHIJJJIJJJJIJIIJIJJIJJIJJIIIJJJJJJJJJJJGHIGGHIIIJHHGEHFFFDEFFECFCEEEDDD PG:Z:MarkDuplicates NH:i:1 HI:i:1 nM:i:0 AS:i:176
FCC0WRYACXX:1:1106:15737:96571 83 X 68206628 255 46M3479N43M = 68201680 -8516 TTCTCATACTGGCAATAGTATTTCACCCAGGATATTCCTAAAGCCCATTTCTCTTGTGTATAGAGATAGCCTTCAATAGTTGGCTGTCC CCCDDCCCDDDCCEDDEDCDEDB9GFFFFGGDGGIIGGEHFAAGFIGGCGDGFCHGIEGGIIHDBEEGIGIJGIIGIIJJGHIHGEFGH PG:Z:MarkDuplicates NH:i:1 HI:i:1 nM:i:0 AS:i:179
But I want only these reads if the fourth column (start position of the reads) in the region. It is possible to get a solution with samtools ?
in my example the start position of the read FCC0WRYACXX:1:1103:14931:24492 is not in the region X:68206675-68208264.
data: RNA-seq paired-end
2 answers
Would this work, without resorting on additional third party tools?
samtools view aln.bam X:100-1000 | awk -v FS='\t' '$4 >= 100 && $4 <= 1000'
If you want also the bam header:
samtools view -h aln.bam X:100-1000 | awk -v FS='\t' '$1 ~ "^@" || ($4 >= 100 && $4 <= 1000)'
(Not tested)
Good call. I'm using a variation of this command in order to pull out SNP genotype-specific reads in order to divide my BAM by region and genotype
How about running bedtools - intersect on your bam file. That should you give you what you want. I just noticed, there is already an example of exactly what you are after on that page.
It says you can use bam as input, but you can't tell bedtools to output a read if the start is inside a target interval. If a read starts before the interval and ends somewhere inside or after the interval, bedtools will print it. The OP wants to discard such reads instead. Did I get it right?
it does not matter which end is obtained by the given list of intervals
Exactly, so if the start of an interval (a read) is outside the interval, the OP wants to reject it. The OP wants to keep reads whose start is inside the interval. So if the read has start and end coordinates chr1:5-20 and the target interval is chr1:10-20, the read has to be rejected, but bedtools will instead return it.
By the way, I'm not considering paired-end reads here.
Log in to answer this question.
-
featureCounts fragment count
written by xosugarxspicexo •
-
Strand specific quantification with strand-splitted BAM files
written by compuTE
-
HTSeq error: local variable 'NH' referenced before assignment
written by senowinski •
-
Samtools view behaviour with skipped region from the reference (N in Cigar)
written by Bastien Hervé
-
Picard Mark Duplicates
written by Varun Gupta
-
Filter bam files using a bed file: Why is the mate missing?
written by komal.rathi
-
Picard Markduplicates Problem
written by xd_d
-
Does htseq-count recognize soft clipping generated by STAR?
written by Dejian
-
How do I configure the read name output from STAR?
written by jgbradley1 •
No, it is not possible with samtools. You will need to write a script and work with the sam file for your purpose.
Given the enormous amount of tools available it's very likely that something exists for this quite specific problem, without the need of writing a custom script. (Although I agree it could be a solution - it's for sure not the easiest for everyone).