Hi,
I have a list of Refseq annotated gene coordinates and I would like to extract the reads for the same in one go. I have read many posts here and at Seqanswers, but those are restricted to a particular region. For example
samtools view input.bam "Chr10:18000-45500" > output.bam
How can one do this for a list of genes with the Transcription start and end positions?
I tried doing this by extracting theit extr Refseq gene names and the associated coordinates from the genome browser at UCSC. Then using samtools view input.BAM "mygenes.bed" > output.bam. But it extracts the reads and the format is not BAM.
I am new to NGS can someone help.
•
written
by
amitpande74 •
3 answers
samtools view -b -L regions.bed input.bam > selected.bam
You should be careful with 0-based and 1-based coordinates
•
written
by
GouthamAtla
Given a sorted file called genes.bed with your Tx starts and ends, you can use BEDOPS bedops and bam2bed running in a bash environment:
$ bedops --element-of 1 <(bam2bed < input.bam) genes.bed > answer.bed
Fast and memory efficient, and useful if you want BED as output for downstream ops.
•
written
by
Alex Reynolds
Log in to answer this question.
-
Filter out Bam not overlapping with Bed File (keeping a read and its mate)?
written by Eliveri •
-
Extract paired-end reads with one end falls within a given region from a bam file
written by Wang Cong •
-
Extract reads overlapping a specific region in bam file
written by User000
-
Converting samtools view <region> output (multiple reads) into single fasta sequence (consensus?)
written by a1ultima
-
What does samflag 1028 mean?
written by junnna •
-
Replace @RG group by the one of another BAM
written by da.de •
-
Extracting reads from given region
written by KVC_bioinfo
-
Samtools gives "truncated file" when trying to view a "reheaded" file
written by olavur
-
How to quickly count the number of reads at a specific coordinate in a bam file
written by klin •
-
Fail to determine the sequence name
written by zwang10 •
