Any tools for extracting aligned part of genome from sam or bam file
1 answer
I can think of a dirty way to achieve what you want. It uses my Sam::Parser from proovread (pacbio correction software). The code below is a draft. It prints a FASTA record corresponding to the aligned part of each read. It cannot handle reverse-complement mappings and probably has some other flaws. But it's an idea. Let me know, if you want to give it a try and I could help you to further improve the snippet.
git clone https://github.com/BioInf-Wuerzburg/proovread.git
export PERL5LIB=proovread/lib/:$PERL5LIB;
BAM=genome.bam
REF=genome.fa
samtools faidx $REF;
samtools view $BAM |
perl -MSam::Parser -e '
my $sp=Sam::Parser->new();
while(my $aln=$sp->next_aln){
next if $aln->is_unmapped();
print $aln->rname,":",$aln->pos,"-",(length $aln->seq_aligned),"\n"
}' |
xargs samtools faidx $REF > ref-reads.fa
Sure, I'm a perl guy, as you can see :). Most of it, you could just do on command line. The only ugly thing is, that you need to parse the cigar string in order to get the reads real length and to account for clipped stuff. I'd be interested to see your solution.
Log in to answer this question.
More posts like this
-
Mapping paired end reads with ngm and samtools, using prefixes and suffixes for creating vcf event…
written by Human •
-
Count expression of 10x scRNA-seq after UMI extraction
written by tomas4482 •
-
Modifying fastq base at specific reference location on different length reads
written by yryan •
-
Sub-sampling RNA-seq data
written by mohammedtoufiq91
-
Match reads from sam (first col) to fastq file and extract to a file
written by sofie_carolina •
-
Remove reads from fastq file
written by varsha619
-
workflow for variant calling using GATK in multiple paired end reads
written by reza
-
Problem using rsem function for validate sam file
written by pier-luc.dudemaine •
-
any tools for updating reference sequence file to include the variants from the vcf file?
written by mheyda •
-
tools for sequence separation
written by mheyda •
Why exactly do you want to do this? If you want to do assembly, then why do you only care about aligned regions? Similarly, why would you do variant calling before-hand?
I don't want to do assembly at this moment, I want to see errors in reads. actually I am working on error correction method for reads, so I want to have the correct read in advance to see how can I correct them.
Ah, your tags were somewhat poorly chosen then. I suspect you'll need to write a tool to do this. I suspect you'll end up needing to rewrite part of the htslib mpileup functionality, simply to allow handling alignments from the stack prior to removing them.