I'm tying to extract a specific region of a bam-file into a fasta-file (ultimately). All of the methods I've tried so far give me all reads that OVERLAP the desired region, I'm trying to find a way to trim those to only the desired region.
I've tried:
samtools view
samtools view compiled.sorted.bam ConB:2185-2195
intersectBed
intersectBed -b test.bed -abam compiled.sorted.bam -ubam > out.bam
but these will give the entire read that overlaps my desired region, I'm trying to get something that will trim everything to sam/bam file where the 'reads' are 10 nucleotides long. Am I just missing a flag somewhere to limit the returned region?
3 answers
I wrote SAM4WebLogo for Sequence Logo For Different Alleles Or Generated From Sam/Bam and I think it could do what you need
$ java -jar dist/sam4weblogo.jar -r seq1:80-110 sorted.bam 2> /dev/null | head -n 50
>B7_593:4:106:316:452/1
TGTTG--------------------------
>B7_593:4:106:316:452a/1
TGTTG--------------------------
>B7_593:4:106:316:452b/1
TGTTG--------------------------
>B7_589:8:113:968:19/2
TGGGG--------------------------
>B7_589:8:113:968:19a/2
TGGGG--------------------------
>B7_589:8:113:968:19b/2
TGGGG--------------------------
>EAS54_65:3:321:311:983/1
TGTGGG-------------------------
>EAS54_65:3:321:311:983a/1
TGTGGG-------------------------
>EAS54_65:3:321:311:983b/1
TGTGGG-------------------------
>B7_591:6:155:12:674/2
TGTGGGGG-----------------------
>B7_591:6:155:12:674a/2
TGTGGGGG-----------------------
>B7_591:6:155:12:674b/2
TGTGGGGG-----------------------
>EAS219_FC30151:7:51:1429:1043/2
TGTGGGGGGCGCCG-----------------
>EAS219_FC30151:7:51:1429:1043a/2
TGTGGGGGGCGCCG-----------------
>EAS219_FC30151:7:51:1429:1043b/2
TGTGGGGGGCGCCG-----------------
>B7_591:5:42:540:501/1
TGTGGGGGCCGCAGTG---------------
>EAS192_3:5:223:142:410/1
TGGGGGGGGCGCAGT----------------
>B7_591:5:42:540:501a/1
TGTGGGGGCCGCAGTG---------------
>EAS192_3:5:223:142:410a/1
TGGGGGGGGCGCAGT----------------
>B7_591:5:42:540:501b/1
TGTGGGGGCCGCAGTG---------------
>EAS192_3:5:223:142:410b/1
TGGGGGGGGCGCAGT----------------
•
written
by
Pierre Lindenbaum
using SamJS
java -jar dist/samjs.jar -e '!( !record.getReadUnmappedFlag() && record.getReferenceName().equals("ConB") && record.getAlignmentStart()>=2185 && record.getAlignmentEnd()<=2195)' input.bam
ha I see, I misunderstood your question. I'll delete my answer.
•
written
by
Pierre Lindenbaum
(1) Read the "README.md" file. (2) This isn't an answer, I'll be deleting it now.
Log in to answer this question.
-
Piping a bam file to samtools calmd
written by rubic
-
Gene variation study with number of genomes
written by Eva_Maria
-
Piping samtools commands
written by jyu429
-
Normalized coverage histograms for selected chromosome region
written by biotech
-
Getting list of chromosome names from indexed BAM file
written by Alex Reynolds
-
Reduce Bam File With Specific Loci In Bedfile
written by marcodpc •
-
Change Notation Of Chromosome In Bam File
written by hellbio
-
Extract Only Paired-End Reads That Map A Specific Interval
written by Nicolas Rosewick
-
Failed Due To Truncated File Or Corrupt Bam Index File Error When Downloading Bam File From Samtools
written by Pascal
-
Counting Number Of Bam Reads Directly Within Set Of Intervals With Bedtools
written by User 9996
So if read spans the boundaries, you want to retrieve just that part of the read that is inside that region?
Correct. I'd prefer to exclude things that are only partially inside the region ... although I can parse that out in my downstream analysis.
There's nothing in samtools, at least, to trim reads at a given boundary, since that's not exactly a common need. I suspect you'll need to code this up yourself (I can foresee some annoyances there).
Yeah, that's what I'm seeing. I figured it would be a more common request, but I guess not ... python to the rescue!