This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to find start and end position of sequence from consensus fasta file?

Hello

I have performed consensus function in samtools and got the consensus fasta sequence from BAM file, and missing positions were replaced with 'N'.

I would like to find my target sequence's start and end position. My sequence is like

>chr1
NNNNAGTATANNNTATGNNNNN

all the heads and ends are N, and my target sequence is in the middle, while some N could be found in the target sequence, I only want to find the start and end of the position where the sequence is not 'N'.

For example, in my case, the start and end position is chr1:5-17 (including three N in the middle of the target sequence).

I have about 1,000 sequences like my example in the same fasta file, Does anyone know how to find the start and end position for each sequence?

Thanks a lot

fasta samtools
• 1,897 views • link
written by 沛煒 • 沛煒

Maybe map the sequences with some aligner like gmap, and then convert the output bam file to bed so that you get the mapping coordinates of your sequence?

• 0 views • link
written by iraun iraun

1 answer

It's assumed here that your fasta file is linear (no linebreaks in sequences). Ugly, but should work..

cat file
>chr1
NNNNAGTATANNNTATGNNNNN
>chr2
NNNNNNNNAGTATANNNTATGNNNNN

paste -d $'\t' \
    <(paste -d $'\t' - - < file) \
    <(paste -d $'\t' \
        <(paste -d $'\t' - - < file | awk 'BEGIN{FS="\t"}{print $2}' | awk 'BEGIN{FS="A|G|C|T"}{print length($1)}') \
        <(paste -d $'\t' - - < file | awk 'BEGIN{FS="\t"}{print $2}' | rev | awk 'BEGIN{FS="A|G|C|T"}{print length($1)}')) \
    | awk 'BEGIN{FS="\t"}{print substr($1,2)":"$3+1"-"length($2)-$4}'
chr1:5-17
chr2:9-21
• 0 views • link
written by 5heikki 5heikki

That works! Thanks for your help

• 0 views • link
written by 沛煒 • 沛煒

Log in to answer this question.

More posts like this