How to merge exons based on gene ids?

I have multiple exons with header gene_id, I wanted to merge the multiple exons sequence to a single string based on header gene_id. The file I have,

   >"gene_id:XLOC_000656"::KB317696.1:3536054-3536229
   AGACCATTTAAATTCTAACTACCAACTGTCACATATAACTCAGCTTCACATAATGTTTGGATCTCGTCGCTCTCGTGCTCCCGCTGCTGGAATGGGTGGCGGAGGTTATGGAGCTGCCCCCCCTCGACGAAGAGGATTTGGAATGGGAGGCGGGCGTCGCGGCGCTGGCTGTGGG
   >"gene_id:XLOC_000656"::KB317696.1:3536284-3536346
   CGCAACCCATGGGAGGGGGAATGGGTGCTGGAAGAACCGGTGGCATGGGGATGGCTCCTGGT
   >"gene_id:XLOC_000656"::KB317696.1:3536407-3536551
   GGGTTGTTGCTCCTCGCGTTAGGACCCGTGATCGCATCAAGCGACTCTTTGGATTTGGCCGTTCGCGTCGTAGGACGGCTATGTACTAAACCAGCCCGTGGTGCAAGCGTGCTTTGAAAGGAAGAAGTGAATAATGTATCATCC
   >"gene_id:XLOC_001528"::KB317696.1:3526472-3526569
   GCTGATCAGAAGTTTGCAGATGATAAAACACGCGCTGCAAAGCTAGAGCAACGCTCTCGGAAGGTATGTCATCCTCGGTAGGTATCTGGTACACGGC
   >"gene_id:XLOC_001528"::KB317696.1:3526622-3526730
   CTTCGGAAATACCGTGTGCAATAGAGCGATTGTAGTAGTGAATTCATATAACATGTTGCACCTTGGTTCTTCAATCCAACGAAACCCGTCTCCTTCTTTGAGTCATAA
   >"gene_id:XLOC_001528"::KB317696.1:3526786-3526969
   TTACAAAGTTATGCCATAGTACACCAGTTGGATCTTTCAGTACTCGAACGAACACGGTTATGTCAGCCTCGTCGCCTTCAATGATCGGTCTACGTCCTTCAATGCCTGTAAAGAGCTTTCGGAGGTCATGGAATCGTGTGAAACCCCAATCGCATTCCTCGGCAATGAACCGATGGTGCGCGT
   >"gene_id:XLOC_001528"::KB317696.1:3527023-3527170
   ACTACACACAAAGTTCGTAGGATCTTGCGGGTTAGACATGACAAGCGCAAATTGAGCGCATGCGTGCCATCCCTCGGGTGCGTTCTTGGGATTCGCATAGTCAAGATATATGGACACTGTATCGTTTGGCGGCGAATTGGAGTTCCC

expected output

  >"gene_id:XLOC_000656"  
  AGACCATTTAAATTCTAACTACCAACTGTCACATATAACTCAGCTTCACATAATGTTTGGATCTCGTCGCTCTCGTGCTCCCGCTGCTGGAATGGGTGGCGGAGGTTATGGAGCTGCCCCCCCTCGACGAAGAGGATTTGGAATGGGAGGCGGGCGTCGCGGCGCTGGCTGTGGGCGCAACCCATGGGAGGGGGAATGGGTGCTGGAAGAACCGGTGGCATGGGGATGGCTCCTGGTGGGTTGTTGCTCCTCGCGTTAGGACCCGTGATCGCATCAAGCGACTCTTTGGATTTGGCCGTTCGCGTCGTAGGACGGCTATGTACTAAACCAGCCCGTGGTGCAAGCGTGCTTTGAAAGGAAGAAGTGAATAATGTATCATCC
  >"gene_id:XLOC_001528"

GCTGATCAGAAGTTTGCAGATGATAAAACACGCGCTGCAAAGCTAGAGCAACGCTCTCGGAAGGTATGTCATCCTCGGTAGGTATCTGGTACACGGCCTTCGGAAATACCGTGTGCAATAGAGCGATTGTAGTAGTGAATTCATATAACATGTTGCACCTTGGTTCTTCAATCCAACGAAACCCGTCTCCTTCTTTGAGTCATAATTACAAAGTTATGCCATAGTACACCAGTTGGATCTTTCAGTACTCGAACGAACACGGTTATGTCAGCCTCGTCGCCTTCAATGATCGGTCTACGTCCTTCAATGCCTGTAAAGAGCTTTCGGAGGTCATGGAATCGTGTGAAACCCCAATCGCATTCCTCGGCAATGAACCGATGGTGCGCGTACTACACACAAAGTTCGTAGGATCTTGCGGGTTAGACATGACAAGCGCAAATTGAGCGCATGCGTGCCATCCCTCGGGTGCGTTCTTGGGATTCGCATAGTCAAGATATATGGACACTGTATCGTTTGGCGGCGAATTGGAGTTCCC

python biopython

• 1,442 views • link

updated 4.8 years ago by Ak •

• written 6.0 years ago by akashbala0 •

Please use the formatting bar (especially the code option) to present your post better.
code_formatting

Thank you!

• 0 views • link

written 6.0 years ago by GenoMax

1 answer

Since this is still unanswered, I'll post an R solution until someone is kind enough to answer with the requested biopython solution.

Read the fasta file in.

library("Biostrings")
library("tidyverse")

fasta <- readDNAStringSet("test.fasta")

You'll end up with a Biostrings object of the fasta file.

> fasta
DNAStringSet object of length 7:
    width seq                                               names               
[1]   175 AGACCATTTAAATTCTAACTACC...GGCGTCGCGGCGCTGGCTGTGGG "gene_id:XLOC_000...
[2]    62 CGCAACCCATGGGAGGGGGAATG...GTGGCATGGGGATGGCTCCTGGT "gene_id:XLOC_000...
[3]   144 GGGTTGTTGCTCCTCGCGTTAGG...AAGAAGTGAATAATGTATCATCC "gene_id:XLOC_000...
[4]    97 GCTGATCAGAAGTTTGCAGATGA...TCGGTAGGTATCTGGTACACGGC "gene_id:XLOC_001...
[5]   108 CTTCGGAAATACCGTGTGCAATA...CCGTCTCCTTCTTTGAGTCATAA "gene_id:XLOC_001...
[6]   183 TTACAAAGTTATGCCATAGTACA...GGCAATGAACCGATGGTGCGCGT "gene_id:XLOC_001...
[7]   147 ACTACACACAAAGTTCGTAGGAT...TTTGGCGGCGAATTGGAGTTCCC "gene_id:XLOC_001...

You can then concatenate the sequences with the same gene id.

new_fasta <- fasta %>%
  {data.frame(seq=as.character(.), name=names(.))} %>%
  mutate(name=str_extract(name, "(?<=\")gene_id:[[[:alnum:]]_]+")) %>%
  group_by(name) %>%
  summarize(seq=str_c(seq, collapse="")) %>%
  {set_names(.$seq, .$name)} %>%
  DNAStringSet

And now you should have a Biostrings object of the new sequences, where the sequences for each gene_id were concatenated based on the order they appeared in the fasta file.

> new_fasta
DNAStringSet object of length 2:
    width seq                                               names               
[1]   381 AGACCATTTAAATTCTAACTACC...AAGAAGTGAATAATGTATCATCC gene_id:XLOC_000656
[2]   535 GCTGATCAGAAGTTTGCAGATGA...TTTGGCGGCGAATTGGAGTTCCC gene_id:XLOC_001528

You can then save this back to a fasta.

writeXStringSet(new_fasta, "new.fasta")

• 0 views • link

written 6.0 years ago by rpolicastro