I have been trying to troubleshoot this for the past couple days but I can't seem to find the issue, any help you could provide is greatly appreciated. I am a noob but will be writing similar scripts to this in the future- so any added resources are welcomed.
I have a large file with a bunch of FASTA sequences and I am trying to extract specific ones based on a list of IDs. While I would like the output to be a collection of whatever ID == FASTA sequences it won't give me an output with results. I know that if I change t to a definite string (ie. 'NM_xxxxx') then it will give a desired output so I suspect the issue is in reading from the csv file.
import re
import sys
import csv
from Bio import SeqIO
output = open("matches.fasta", "w")
x = SeqIO.parse(open(sys.argv[1], 'rU'),'fasta')
y = csv.reader(open(sys.argv[2], 'rU'))
for row in y:
print str(row[1])
t = str(row[1])
for fasta in x:
name, sequence = fasta.id, fasta.seq.tostring()
searchObj = re.search( t, fasta.id, re.M|re.I)
if searchObj:
SeqIO.write(fasta, output, "fasta")
output.close()
6 answers
Try pyfaidx, it is really easy to use.
pip install pyfaidx
then to extract your seq
xargs faidx -d "<here the separator used in the IDs of your fasta>" <Your fastafile.fasta> < <Here the file containing the IDs>
Example
xargs faidx -d "|" My_seq.fasta < My_Ids.txt
PS: this is a clear example of the great "teaching power" of Biostars. extract sequences from multi fasta using partial ID
What you are missing is seek() in biopython. Assuming the seq_ids in a file called "names"
How to use seek() with Biopython SeqIO object
from Bio import SeqIO
fasta=SeqIO.parse("fasta.fa","fasta")
seq_dict={}
for record in fasta:
seq_dict[record.id]=record.seq
# If you feel its memory intensive, try next solution
for line in open("names","r"):
line=line.strip()
if line in seq_dict.keys():
print ">"+line+"\n"+seq_dict[line]
Otherwise write a small function like:
from Bio import SeqIO
def getSeqRec(seq_id):
fasta=SeqIO.parse("fasta.fa","fasta")
for record in fasta:
if seq_id==record.id:
return record
for line in open("names","r"):
line=line.strip()
seq_rec=getSeqRec(line)
print ">"+seq_rec.id
print seq_rec.seq
Does it need to be Python? samtools faidx does it.
samtools faidx Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa chr22
•
written
by
mark.ziemann
In linux
#exract fasta sequences with a reference list ALL SEQUENCE LENGTH
#LIST is your list as text and IN is your fasta
awk 'BEGIN{while((getline<"LIST")>0)l[">"$1]=1}/^>/{f=l[$1]?1:0}f' IN.fa > OUT.fa
Here is a script I wrote. Before using it, Perl and perl module BioUtil should be install.
fasta_extract_by_pattern.pl could extract FASTA sequences by header or sequence, exactly matching or regular expression matching are both supported. The query pattern could read from files. And negation of the result is also easy to get. What's the most important, it could read from STDIN.
Combining fasta2tab and tab2fasta with cvs_grep could also have the same function.
Examples
1. sequences WITH "bacteria" in header
fasta_extract_by_pattern.pl -r -p Bacteria *.fa > result.fa
2. sequences WITHOUT "bacteria" in header
fasta_extract_by_pattern.pl -r -n -p Bacteria seq1.fa seq2.fa > result.fa
3. sequences with TTSAA (AgsI digest site) in SEQUENCE. Base S stands for C or G.
fasta_extract_by_pattern.pl -r -s -p 'TT[C|G]AA' seq.fa > result.fa
4. sequences (read from STDIN ) with header that matches any patterns in list file
zcat seq.fa.gz | fasta_extract_by_pattern.pl -pf name_list.txt > result.fa
See more: Manipulation on FASTA/Q format file
•
written
by
shenwei356
To extract sequences with specific names from a large file efficiently, I suggest you use my filterbyname tool:
A: Extract sequence with header from a fasta file with specific ID given in another
Log in to answer this question.
-
Tool: Converting Degree Minutes (DM) to Decimal Degrees (DD) using python script
written by kuttibiotech2009 •
-
How to convert dataframe properly to anndata?
written by bioinfo •
-
How to get rid of all copies of repeated lines in BLAST result text files
written by Hau •
-
uploding VCF file to CADD and Ensemb
written by Eliza •
-
finding overlapping motifs to increase length of motif
written by andrea •
-
10X Genomics, Single Cell Multiome ATAC-seq and RNA-seq data
written by yzhao140 •
-
get organism name from fasta file header (PYTHON)
written by biohacker_tobe •
-
python script to sort blast results and fetch best matches from NCBI
written by skbrimer
-
problem in generating 6 frames
written by swats1994 •
-
Help: do I always have to re-open sequence file within a for loop in python?
written by pawlowac
I hope you read through the dozens of ID matching posts here.
I did, but the fact my experience is almost exclusively python and very limited at that means I am unable to take advantage of related posts.
From your description, there could be one problem with the data that is a bit undetectable unless you look for it. For a few sequences, after reading them with Bio.SeqIO, try printing their
ID,nameanddescriptionattributes. Check if your expectations of what should be printed are in fact what are printed.If the FASTA header has white spaces, the SeqIO parser can be a bit crazy. I'd usually say that the part before spaces goes into id and the part after goes into description, but this assumption was recently proven wrong.
Check out if the data conforms to expectations, and then dig into your code for changes.
The main issue is that the FASTA header has a wierd header (ie. >NM_xxxxxxx_2000_from_chromosome_spongebob) which is why I was trying to use RE to match with my candidate list containing just the RefSeq Ids (NM_xxxxxx). Else, Bio.SeqIO is behaving like it should given the control tests I've run.
Ah, I see. Have you tried printing
str(row[1]), as plain text as well as an ASCII sequence? Maybe a new line or a black space at the end (or a weird unicode character) is messing up your RE search.You were right, looks like the export .txt left some weird ascii characters. Thanks for the help Ram!
In linux
#exract fasta sequences with a reference list ALL SEQUENCE LENGTH
#LIST is your list as text and IN is your fasta
>awk 'BEGIN{while((getline<"LIST")>0)l[">"$1]=1}/^>/{f=l[$1]?1:0}f' IN.fa > OUT.fa
Couple of notes:
a. LIST is a file, not text. substituting that can be a bit non-trivial
b. Heng Li's bioawk is awk for biological formats. It will obviate the need for multi-line checks and header line awareness in your awk script.