1 answer
Have a look at https://github.com/shenwei356/seqkit and its subseq command: "get subsequences by region/gtf/bed, including flanking sequences". All you need is fasta file of you genome of interest and a gff file (gene annotations). Command could look like this:
seqkit subseq --up-stream 200 --downstream 200 --gtf genome.gff genome.fa
Yes, something like the snippet below could work. Needs one little extra piece of command line magic, though. I'm assuming you are working on a bash-like command line:
# get orfs with flanks as before
seqkit subseq --up-stream 200 --downstream 200 --gtf genome.gff genome.fa > all-orfs-with-flanks.ffn
# get orfs w/o flanks
seqkit subseq --gtf genome.gff genome.fa > all-orfs-without-flank.ffn
# grep only sequences that end with TAA and get their ids
seqkit grep -rsp 'TAA$' all-orfs-without-flank.ffn | perl -ne 's/^>// && print' > orfs-with-TAA-stop-ids.txt
# extract those sequences from the orfs with flanks
seqkit grep -f orfs-with-TAA-stop-ids.txt all-orfs-with-flanks.ffn > orfs-with-TAA-stop-and-flank.ffn
Hey, thackl All commands except the the last one worked perfectly well. After playing a bit with the sequences it seems to me that there is a subtle matching problem between the sequence ids on the ffn file and the ones in the list orfs-with-TAA-stop-ids.txt I manually copied/pasted a couple of sequences from the all-orfs-with-flank.ffn file into a new file and the corresponding ids into another file and then the 'seqkit grep -f' command worked out. The funny thing is that the ids in the original all-orfs-with-flank.ffn and orfs-with-TAA-stop-ids.txt files look identical. Any idea how to fix this?
OK, this worked for me on a dummy example. Minor fixes to old code & strip the flank annotations from IDs. Only problem, the final file doesn't contain "gene names" but id as chr:start-end...
seqkit subseq --up-stream 200 --down-stream 200 --gtf genome.gff genome.fna > all-orfs-with-flanks.ffn
seqkit subseq --gtf genome.gff genome.fna > all-orfs-without-flank.ffn
seqkit grep -rsp 'TAA$' all-orfs-without-flank.ffn | perl -ne 's/^>//;s/ $//; print' > orfs-with-TAA-stop-ids.txt
perl -pe 's/>(\S+_\d+-\d+:.).*/>$1/' all-orfs-with-flanks.ffn |
seqkit grep -f orfs-with-TAA-stop-ids.txt > orfs-with-TAA-stop-and-flank.ffn
Log in to answer this question.
-
Predict species from DNA sequences
written by p2016k •
-
Given set of genomic sequences find potentially enriched genes?
written by chipolino •
-
Trim 3' ends from single-end reads
written by Maria •
-
Extract flanking regions from bam files based on coverage threshold
written by laura.bertola •
-
Is there any tool to calculate dn/ds for orthologues?
written by ThulasiS
-
How to retrieve gene sequences from whole genomes
written by nguyen.lamsyd •
-
pull out flanking sequences of snp
written by myprogramming2016 •
-
retrieve orthologous sequences from human fasta file
written by jfertaj
-
Identify common genes and their respective DNA sequences
written by sam •
-
How To Retrieve Upstream Sequences Of A Single Gene From Multiple Bacterial Species?
written by Magdalena •
