I have a BAM file that does not contain any chromosome or position information, our file contains only scaffold of a new genome from specie B. I'm looking to find scaffolds/reads matching a gene (35kb long) from specie A in the BAM file of specie B. I expect some differences (at the SNP level). Since our bam file doesn't contain chromosome and position information, is there a way to tell samtools, this my sequence, find any reads/scaffold that can match (with at least 60% of similarity) that sequence?
PS: Thanks for answering thou I do not have access to BWA so far.. otherwise I would have done it with BWA long time ago.. I work remotely on a server with limited capability T.T
OUR BAM FILE THAT CONTAINS ONLY SCAFFOLD INFO (NO CHROMOSOME OR POSITION)
FCC00CUABXX:7:2105:3995:105611#NNNNNCGT 73 scaffold12 24440 37 100M = 24440 0 AACCATTTCGTTTCTGAATGAAGAATACAATACTGGAAATGAAAACTTCACTACAGGAACTCAATAGCAAAGTAGATGATACAGAAGAACAGATCAGCAA HHHHHHHHHHHHHHHHHHHHHHHHFHHHHHHHHHHHFHHHHHHHHHHHHHHHHHHHFHFHHHHHHHHHDFBEBFFHFHHHHH9HDFFFEHHHHHDDFE<E X0:i:1 X1:i:0 MD:Z:100 RG:Z:DONdtxDABDMAAPE XG:i:0 AM:i:0 NM:i:0 SM:i:37 XM:i:0 XO:i:0 XT:A:U
FCB01AKABXX:1:1202:4275:93484#TCNGNGCN 97 scaffold12 24443 37 49M = 21972 -2422 CATTTCGTTTCTGAATGAAGAATACAATACTGGAAATGAAAACTTCACT HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH X0:i:1 X1:i:0 MD:Z:49 RG:Z:DONdtxDAADWAAPE.1 XG:i:0 AM:i:37 NM:i:0 SM:i:37 XM:i:0 XO:i:0 XT:A:U
3 answers
convert the reads back to FASTQ and map them with bwa to your gene ?
samtools does not align sequences (as Pierre alludes to); as the name implies, it is used to analyze and process SAM/BAM data, which is assumed to be aligned already. You need to either extract the FASTQ or (if the aligner accepts BAM as input) use it directly and align to sequence B.
EDIT: by 'as Pierre alludes to', I meant his indication that bwa (an aligner) would be required.
•
written
by
Chris Fields
Agreed, see my clarification.
If you could extract your sequences to a fasta file, I have a tool that would do exactly what you are looking for
It uses BLAST, a reference(s) and query sequences, your fasta of extracted sequences. You may specify both identity and alignment length. In the output you'll have the extracted genes of interest in a fasta file
Log in to answer this question.
-
How "bases mapped (cigar)" value in samtools stats is calculated?
written by Tian •
-
Filling gaps in BAM file
written by hemr3 •
-
extract discordant,split reads from BAM file
written by 9521ljh •
-
Get gene counts from *.sam file and annotation
written by tim.ivanov.92 •
-
Trimmed NGS exome data failed to be mapped with bwa
written by yanqiangchang •
-
Can't convert paired end BAM to bed using bedtools
written by colin.kern
-
Modifying bam file sequence names according to 2nd column
written by jyu429
-
BWA: improper orientation but listed as properly paired
written by richard.deborja
-
Should We Dump Illumina Pair-End Mapping Results In Sam With Mapq=0, But Good Template Length
written by PeterPan •
-
Bam File Output Format .. Don'T Understand!
written by madkitty
