1
$\begingroup$

I have a fasta file (gene.fasta)with sequences with names patterns:

>Hsc_gene_9627.t1
ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACAACATCTGCTGTGCAGCCGCTTCA
TTGCTCATTTTGAATGCCGTTACATTGGAGAAGGATGCTAATGATTATGCCGTTGGCGAT

and I have the Ids in text file (seq.txt) which are not the exact as in the fasta file:

HSC_gene_996
HSC_gene_9734
and some of the names came as 
HSC_gene_996|HSC_gene_9734 

How can extract the sequences? I am not experienced in python so please use python for dummies language :) Thank you in advance

$\endgroup$
7
  • $\begingroup$ Does it have to be python? There are far easier ways using only bash, etc. $\endgroup$
    – NatWH
    Commented Jun 3, 2020 at 11:47
  • $\begingroup$ I just have python in my computer, if there is please help me, I tried to find some online tools but it didn't work $\endgroup$ Commented Jun 3, 2020 at 12:37
  • $\begingroup$ You may want to try Joe Healey's bioinfo-tools repository. This script looks like the one to use. $\endgroup$
    – Ram RS
    Commented Jun 3, 2020 at 14:07
  • $\begingroup$ To second @RamRS, that script should work indeed - provided your IDs are exact full or substring matches (excluding the >). It's been a while since I used it, so if it could be tweaked to be more useful, I'm happy to take suggestions. Whilst I love bash etc., I'd probably disagree with NatWH that it is easier than using something like biopython which is purpose built for this. $\endgroup$
    – Joe Healey
    Commented Jun 3, 2020 at 15:08
  • $\begingroup$ OP's problem seems to be their platform. Python excels in being a cross-platform solution, hence OP looking for Python based solutions makes sense. $\endgroup$
    – Ram RS
    Commented Jun 3, 2020 at 15:24

1 Answer 1

1
$\begingroup$

For example, the fasta are like these:

gene.fasta:

>Hsc_gene_9627.t1
ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACA
ACATCTGCTGTGCAGCCGCTTCA
>Hsc_gene_9627.t2
GGGGTTTTCCCC
>Hsc_gene_962.t1
AGTCGTCAGTCAGTAGTCGC

seq.txt:

HSC_gene_9627
HSC_gene_9734 
HSC_gene_996|HSC_gene_9734

We use a few modules, read in the fasta:

import re
from Bio import SeqIO                                                                 
records = list(SeqIO.parse("gene.fasta","fasta"))

Read in the list:

seq = [i.rstrip().lower() for i in open('seq.txt').readlines()]

Here I use two functions, one to write fasta files, and the other to format the header for comparison:

def fasta_out(rec):
    return(">"+str(rec.description)+"\n"+str(rec.seq))

def format(rec):
    return re.sub("[.][^ ]*$","",rec.description).lower()

Define the records to keep:

keep = [fasta_out(i) for i in records if format(i) in seq]

Write them out:

f=open('matched.fasta','w')
f.writelines("\n".join(keep))
f.close()

We can check:

open('matched.fasta').readlines()                                                          
Out[50]: 
['>Hsc_gene_9627.t1\n',
 'ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACAACATCTGCTGTGCAGCCGCTTCA\n',
 '>Hsc_gene_9627.t2\n',
 'GGGGTTTTCCCC']

Most likely there are some simple while, for loops that does not require modules..

$\endgroup$
1
  • $\begingroup$ Thank you very much @StupidWolf $\endgroup$ Commented Jun 10, 2020 at 12:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.