Thanks. But, I want to combine with pyfasta so that I can extract gene from each sample file + add sample name (file name) to header in one step. Can you suggest some easier way, such as awk, sed etc. As awk is working for one file and I did not do any linerization eventhough that file is multi-line fasta.
Add filename to fasta headers of multiple fasta files inside loop
Hi,
I have 10 fasta files (each file with 20 gene sequences from each of the 10 samples). I would like to create 20 files, specific to each gene from 10 samples. I proceeded as follows to extract genes with the file_name in header:
pyfasta extract --header --fasta test.fasta gene_name1 | awk '/^>/ {$0=$0 "_sample1"}1' > gene_name1.fasta
Output:
>gene_name1_sample1
ATGC
I am successful in creating multiple gene fasta files for each gene from each sample (a part from loop):
pyfasta extract --header --fasta $sample.fasta gene_name1 >> gene_name1.fasta
pyfasta extract --header --fasta $sample.fasta gene_name2 >> gene_name2.fasta
But, I am unable to add file_name to the header of files in loop (but can do for 1 file as mentioned in the beginning).
Kindly guide.
Thanks.
2 answers
The question is pretty generalized with the use of terms. From what I understand, this should work:
If each sample file are not single line fasta, linearize, on commandline:
for file in *.fasta; do awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' < $file > "`basename $file .fasta`_single-line.fasta"; done
Then run this python script:
#!/usr/bin/env python
import glob
from collections import defaultdict
genes = defaultdict(list)
for file in glob.glob('*_single-line.fasta'):
with open(file, 'r') as f:
for line in f:
if line.startswith(">"):
# Gene ID : (sample name, sequence)
genes[line.strip().split('>')[1]].append((file.split('.fasta')[0], next(f).strip()))
for g in genes:
# open fasta per gene name
with open(g + '.fasta', 'w') as out:
# for each sample name and sequence
for n in range(len(genes[g])):
# print to file, >gene_sample, sequence
out.write('>' + g + '_' + genes[g][n][0] + '\n' + genes[g][n][1])
I want to combine with pyfasta so that I can extract gene from each sample file + add sample name (file name) to header in one step
OK, so this is clearly a two step process, however the code can be converted to not have to use the awk statement first, and simply input the multi-line fasta.
Let's first check if the above does what you need.
Put this into a bash script, run.sh:
#!/usr/bin/bash
for file in *.fasta; do
awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' < $file > "`basename $file .fasta`_single-line.fasta"
done
python combine.py
And the above python code into, combine.py.
Run the pipeline as bash combine.py in the directory with your sample fastas. This will run in one step.
Thanks @st.ph.n. run.sh worked. :)
combine.py also worked but file name did not appear near to 1st sequence and no sequence for last. Also, all files are of same size, but I did not understand that all files are in order except file 8 and 9.
#gene1.fasta
>gene1_file10
ATGC...............................>gene1_file8
ATGC...............................>gene1_file9
NTGC...............................>gene1_file7
ATGC...............................>gene1_file6
NNGC...............................>gene1_file5
ATGC...............................>gene1_file4
ATGC...............................>gene1_file3
ATGC...............................>gene1_file2
ATGC...............................>gene1_file1
ATGC...............................
#gene2.fasta
Same for gene2 and other genes
I don't see why the file is appearing like you posted. I ran a test on this script, and it appears fine. What python version are you using? python --version
try changing this line:
out.write('>' + g + '_' + genes[g][n][0] + '\n' + genes[g][n][1])
to this, to print to stdout and file:
print '>' + g + '_" + genes[g][n][0], '\n', genes[g][n][1]
print >> out, '>' + g + '_' + genes[g][n][0], '\n', genes[g][n][1]
You just need to add a '\n' at the end of the string using the write method.
• 0 views
•
link
written
by
Matt Shirley
Thanks Matt Shirley
[int] specifies in index, in this case line.strip().split('>')[1] is splitting the header into '>' and 'gene name', and python starts at 0 so we're taking the gene name. Sae with splitting the file name into parts, 'file name' and '.fasta', and we're taking the file name.
You can use pyfaidx to do this in one step:
$ pip install pyfaidx
$ faidx -x -e "lambda x: x + '_sample1'" test.fasta
The faidx script is comparable to the pyfasta script, but has many more features. It also creates/uses the .fai index, just like samtools. The -x argument splits the input file into separate files, named using the header. The -e argument takes a python lambda function that returns a modified header name, allowing you to do complex transformations. If you don't pass any region positional arguments, faidx will operate over all sequences in the input FASTA file, so you don't need a bash loop.
• 0 views
•
link
written
by
Matt Shirley
@Matt Shirley Thanks, it will take only one fasta file at a time. My aim is to extract the genes with similar gene name from all the fasta files and make gene specific fasta files with updated header including gene name and file name (so that I should know from which file that gene came) + append the gene sequences in the file with that gene name. Here are the sample input and output files:
Input files:
#file1.fasta
>gene1
ATGC
>gene2
ATGA
>gene3
ATGTTT
#file2.fasta
>gene1
ATGG
>gene2
ATGC
>gene3
ATGTT
Expected output files:
#gene1.fasta
>gene1_file1
ATGC
>gene1_file2
ATGG
#gene2.fasta
>gene2_file1
ATGA
>gene2_file2
ATGC
Log in to answer this question.
More posts like this
-
Add gene expression boxplot from plotCounts to each row of oncoprint via ComplexHeatmap
written by bioinfo8
-
Extract sequence from Fasta using header
written by Princy •
-
Extract the sequence
written by zhichusun •
-
How to save nested for loop data in R
written by archana.bioinfo87
-
Labelling sequences within fasta files according to sample name.
written by Mitra •
-
Extract consensus dna sequence of multiple genes from multiple bam files in fasta format
written by bioinfo8
-
Retrieval of similar positions from samples in a vcf but different from reference
written by bioinfo8
-
CoverageView package_Error in 1:nrow(grDF) : argument of length 0
written by bioinfo8
-
error in fasta sequence extraction based on gene ID
written by Bioinfonext
-
Iterating Through Fasta Sequences With Pyfasta Python Module
written by Eric Normandeau
Can you post an example of some of the headers from one of the sample fasta? Is each sample fasta a multi-line or single line fasta? Is each header, per gene, identical across sample fastas?
Header from sample1.fasta file
Like this upto 20 gene names and their sequences. Each sample fasta files are multi-line. Gene name is same across sample fastas.