This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to count a sequence and cut in the other sequence?

I have a reference sequence (CDS) and an aligned sequence in the same file. Format fasta.aln or aln.The alignment was done with MAFFT.

Input:

RefSeq - - - - - - AAGCTGC

Seq1 AAAAAAGGGGGG

Output I would like:

Seq1 GGGGGG

It would be a removal from sequence 1 according to the symbol "-" of RefSeq. I would like to extract only the CDS after the alignment. Is there any way to do this from the command line or some programming language? I tried to do it with biopython but was not successful!

f = open('Denv4-X-gb_AY947539.txt', 'r')
con = f.readlines()
con = [i.strip() for i in con]
length = len(con[0].split("-")[0])
result = f'{con[0].split("-")[0]} {con[0].split("-")[0][length:]}'
print(result)
f.close()
f = open('Denv4cds.txt', 'a') f.write(f'\n{result}')

Do you have a script or module or library that can do this?

I'm using windows wsl.

I'm a beginner in bioinformatics

the generated file just printed the first line ...

alignment sequence gene software error genome
• 550 views • link
written by USER • USER

I think you need to check the formatting of the post here. It's not clear to me from your description how GGGGG relates to the RefSeq sequence?

I interpret this to mean you want to trim off sequence in your query which corresponds to gap characters in the RefSeq aligned string?

• 0 views • link
written by Joe Joe 
> >lcl|NC_001477.1_cds_NP_059433.1_1/1-10179 [gene=POLY] [locus_tag=DV1_gp1] [db_xref=GeneID:5075725] [protein=polyprotein]
> [protein_id=NP_059433.1] [location=95..10273] [gbkey=CDS]
> ------------------------------------------------------------------------
> ----------------------atgaacaaccaacggaaaaagacgggtcgaccgtctttcaatatgctgaa acgcgcgagaaaccgcgtgtcaactgtttcacagttggcgaagagattctcaaaaggattgctttcaggcca
> aggacccatgaaattggtgatggcttttatagcattcctaagatttctagccatacctccaacagcaggaat
> tttggctagatggggctcattcaagaagaatggagcgcttaacagcacgagccacctgggccac
> caacatacaagtggccataaaccaagtgagaaggctcattgggaatgagaattatctagacttcatgacatc
> aatgaagagattcaaaaacgagagtgatcccgaaggggcactctggtaa-----------------------
> ------------------------------------------------------------------------
> ------------------------------------------------------------------------
> ------------------------------------------------------------------------
> 
> 
> >gb:AB189120|Organism:Dengue/1-10735 virus 1|Strain Name:98901518 DHF DV-1|Segment:null|Subtype:1-IV|Host:Human
> agttgttagtctacgtggaccgacaagaacagtttcgaatcggaagcttgcttaacgtagttctaacagttt
> tttattggagagcagatctctgatgaacaaccaacggaaaaagacgggtcgaccgtctttcaatatgctgaa
> acgcgcgagaaaccgcgtgtcaactgtttcacagttggcgaagagattctcaaaaggattgctttcaggtca
> aggacccatgaaactggtgatggcttttatagcattcctaagatttctagccatacctccaacagcaggaat
> tttggctagatggggctcattcaagaagaatggagcgattaaagtgttacggggtttgtaagagctatgctgcct
> gtgagccccgtccaaggacgtaaaatgaagtcaggccgaaagccacggtttgagcaagccgtgctgcctgta
> gctccatcgtggggatgtaaaaacctgggaggctgcaacccatggaagctgtacgcatggggtagcagacta
> gtggttagaggagacccctcccgaaacataacgcagcagcggggcccaacaccaggggaagctgtaccctgg
> tggtaaggactagaggttagaggagaccccccgcataacaataaacagcatattgacgctgggagagaccag
> agatcctgctgtctctacagcatcattccaggcacagaacgccagaaaatggaatggtgctgttgaatcaac
> aggttct-----------------------------------------------------------------
> ------------------------------------------------------------------------

it's something like this.The GGG part relates to the result of the sequence after cutting. I tried to use a fasta editor but there are many sequences
• 0 views • link
written by USER • USER

Please use the code formatting button (looks like 101010). I tried to format your output myself but it's all over the place.

• 0 views • link
written by Joe Joe

I have done the changes

• 0 views • link
written by USER • USER

Hello projetoic!

This question has been asked again. Please do not open multiple threads for the same question. How do I count a number of ("-") in a reference string and remove it in a second string?

For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.

If you disagree please tell us why in a reply below, we'll be happy to talk about it.

Cheers!

• 0 views • link
updated by GenoMax GenoMax • written by Joe Joe

1 answer

You could get this done by using the -nogaps option of trimAl.

• 0 views • link
written by Mensur Dlakic Mensur Dlakic

But I would not like to remove all gaps ... And according to the gaps of the reference sequence (lcl|NC_001477.1_cds_NP_059433.1_1/1-10179) , remove the characters or sequence of the second sequence (gb:AB189120|Organism:Dengue/1-10735 virus 1) of the organism or species. Is this option the most correct?

• 0 views • link
written by USER • USER

Log in to answer this question.

More posts like this