Hello All,

I have aligned the query sequence with the human genome. I am interested in looking at exons of a specific gene. When I look at the IGV I can see several reads that are mapping to the gene of interest across all the exons. I am confused how can I get the read length of each read which is aligning to exon 1, exon 2 and so on.

I tried to extract the reads giving coordinates of exons using samtools. However, when i converted the resulting bam file to fasta and tried to calculate the read length this gives me the original read length and not the read length of read matching to the region I specified.

I used the following command:

For extracting the reads from original bam file for the region of interest:

samtools view -h in.bam "chr3:region1-region2" > exon.bam

Converting bam to fasta :

samtools bam2fq exon.bam | seqtk seq -A > exon.fasta

for counting the lengths of reads from fasta file:

awk '/^>/ {if (seqlen) print seqlen;print;seqlen=0;next} {seqlen+=length($0)}END{print seqlen}' exon.fasta

Am I doing anything wrong here or missing anything totally? Thank you.