i have a file named gene.txt contain sequence coordinates such as :
orf1=456-12512
orf2=2869-6898
another file contain multiple fasta sequence called project.fasta
i want to extract header with sequence for each one in a separate file, i tried the following seqkit command but it gives invalid std in
awk '{print $0}' gene.txt | seqkit grep -r -n -p $orf1 project.fasta
•
written
by
Mohamed •
1 answer
Please use seqkit faidx which accepts a file containing a list of regions in the format of seqid:start-end similar to samtools faidx.
$ cat gene.txt
orf1=5-7
orf2=1-10
$ cat project.fasta
>orf1
AAAAAACCCCCC
>orf2
TTTTTTGGGGGG
$ seqkit faidx project.fasta -l <(sed 's/=/:/g' gene.txt)
[INFO] 2 patterns loaded from file
>orf1:5-7
AAC
>orf2:1-10
TTTTTTGGGG
when i type the content of gene.txt on terminal without streaming the file in the pipeline and then run the seqkit command line it works without changing the header of fasta file, but i need to pipe all content of gene.txt because it contain a lot of files it will be difficult to change header in fasta file to match both files
yes for @shanwei showed this example it works well beacuse both header matched, but because of header in multifasta file isnot the same so it gives [WARN] sequence not found. as illustrated bellow
$ cat gene.txt
orf1=5-7
orf2=1-10
$ cat project.fasta
>NC_00001.1:456-12512
ATAGATAT
>NC_00001.1:2869-6898
TTTTTTGGGG
so i need to pickup each sequence with header in separate file
but because of header in multifasta file isnot the same so it gives [WARN] sequence not found.
Basic requirement for @shenwei's solution is that gene.txt should have fasta headers that fully match what is in your fasta file. Otherwise there is no point in using his solution because it will not work.
so i need to pickup each sequence with header in separate file
I am not 100% certain what you are saying here. If you have more than one fasta file to pick these sequences from then you need to run the solution for each file (unless you want to cat the files into one giant source file).
Then delete the locations before using seqkit faidx.
$ seqkit seq -i --id-regexp '^(.+?):' project.fasta -o tmp.fasta
or
$ sed -E 's/:.+//' project.fasta -o tmp.fasta
Ideally, just add --id-regexp '^(.+?):' to the seqkit faidx command.
seqkit faidx --id-regexp '^(.+?):' project.fasta -l <(sed 's/=/:/g' gene.txt)
I've just fixed the code, the --id-regexp did not work for faidx before. You can compile from the source or directly download the binary file
Log in to answer this question.
-
How to extract headers from protein fasta file for seqkit analysis?
written by marongiu.luigi
-
find and count motif in FASTA file
written by Matteo Ungaro •
-
Redirect Stdout From twoBitToFa in Bash
written by Laura •
-
vcf to SNPmatrix through vcftools --012 produces extra column
written by liyong •
-
Extracting the coordinates for a specific base from multifasta
written by Louis Kok •
-
Forum: Awk in Bioinformatics
written by Shicheng Guo
-
Problem with FASTA filtering
written by Janey •
-
Extract 6mer and headers from fasta and compile into list
written by omansn •
-
Comparing fasta headers to set of 3 letter codes that is found in the header
written by pasanfernando87 •
-
Extraction Of Fasta Sequence By Nucleotide Position
written by rosarylimyt •
Are you looking to pull out entire sequences or just the coordinate boundaries that are in the examples above?
grepwould not be the correct tool for the former.You may want to use
samtools faidxinstead.