Hi,
I have tried a "home-made" method in order to extract the number of nucleotides per read out of a SAM file but is not working correctly for all the reads due to some deletions (I guess).
grep -e "pattern" my.sam | cut -f 2,3,4,5 > output.txt
To explain, I search for reads with a pattern and then extract some information from every read such as Chromosome, Position, Strand, and Sequence. Then I use R software to count the number of characters of the "Sequence" column and I get the number of nucleotides per read. However, the sequence sometimes might contain a deletion "-" which counts as a character and I get some misplaced reads. I don't want to get rid of these reads. My alignment parameters allow only 1 mismatch so I expect to get a single deletion or addition if that matters.
Is there any way to use the cigar strings of each read and extract the correct number of aligned nucleotides per read?
Thanks for your time, Ioannis
•
written
by
ioannis •
2 answers
using bioalcidaejdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html and looping over the cigar string:
java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar()!=null).forEach(R->println(R.getReadName()+"\t"+R.getContig()+"\t"+R.getStart()+"\t"+R.getCigar().getCigarElements().stream().filter(CE->CE.getOperator().isAlignment()).mapToInt(CE->CE.getLength()).sum()));' input.bam
Hi, you can use the raw output from samtools mpileup (without options), and on the 5th column you have the aligned nucleotides at every position.
Get all the nucleotides from a bam file that is aligning to specific position in the genome
the question was "number of aligned bases per reads" not "number of reads per position"
Log in to answer this question.
-
Identifying a true viral genomic deletion, how to distinguish real deletions from quasispecies nois…
written by tomas •
-
Regex in umi-tools doesn't identify pattern in non-beginning positions
written by kalaivani •
-
(STAR) What causes a high number of reads that are unmapped because they are too short?
written by becky.reese •
-
STAR aligner can't map too short reads
written by Assa Yeroslaviz
-
RNA Seq analysis
written by Micro_bioinfo •
-
ATAC-seq problem with STAR - Multimapping reads
written by Rafael Soler •
-
Terrible RNA mapping result by STAR
written by 610225668 •
-
Finding motifs across a genome
written by ioannis •
-
How can I split the soft clipped reads and map the splitted reads again.
written by fatima.m.zare •
-
Long Read Length, yet STAR says many reads too short
written by dec986
that's not a SAM file. The specification doesn't allow a hyphen. https://samtools.github.io/hts-specs/SAMv1.pdf
I see, but maybe I didn't make it clear. In the Sequence column I do not get any "-" or any other symbols but I get some aligned reads with chromosome and position that do not much the genomic sequence by one nucleotide.
Example: Genomic sequence: AGTCTAGTACCC Aligned sequence: AGCTAGTACCC
In that case there is a single deletion and I wonder if within the SAM file this information can be found and somehow extracted.