No need for open(FNAME, "rU") or IFILE.close(). SeqIO handles all opening and closing itself.
You can modify this simply to:
for record in SeqIO.parse(sys.argv[2], "fasta"):
Hi everyone!
I would like to know how I can extract the 1000 last nucleotides from sequence in a FASTA file and their headers using linux command line.
There are about 2000 sequence in this FASTA file and need to extract the last 1000 nucleotides from each sequence keeping their original headers.
I am a beginner at bioinformatics and I am still learning how to manipulate FASTA files, and I'm stuck with this problem and I can not solve it.
Thanks!
This requires Biopython. Save it as extractEnds.py and chmod +x extractEnds.py to make it executable and then just run it
#!/usr/bin/python
import sys
from sys import argv
from Bio import SeqIO
if (len(sys.argv) < 2):
print("Usage: extractEnds.py file.fna")
sys.exit()
else:
FNAME = argv[1]
IFILE = open(FNAME, "rU")
for record in SeqIO.parse(IFILE, "fasta"):
print(">%s\n%s" % (record.description, record.seq[-1000:]))
IFILE.close()
I'm quite a Python noob so if there's a way to make this better or there are mistakes, please let me know..
No need for open(FNAME, "rU") or IFILE.close(). SeqIO handles all opening and closing itself.
You can modify this simply to:
for record in SeqIO.parse(sys.argv[2], "fasta"):
Thanks. What's the difference between (argv[1], ..) and (sys.argv[1], ..)?
edit. nevermind, found it, it just depends on what gets imported in the start..
edit2. For years I sticked with epic GNU coreutil oneliners, but I'm starting to see the advantages of biopython. Especially if you write scripts that someone else inherits some day.. they will be much happier if you used biopython
If you want it to be robust, or you want the code to be used by anyone but yourself (especially if they're a novice), a dedicated parser like BioPython which covers lots and lots of edge cases is worth its weight in gold.
You could write that code in a reasonably concise one-liner too (this one also has the benefit of taking a second CLI arg for any length of sequence):
python3 -c "import sys; from Bio import SeqIO; [print(f'>{r.id}\n{r.seq[-int(sys.argv[2]):]}') for r in SeqIO.parse(sys.argv[1],'fasta')] ;" seqs.fa 1000
(Python 3 only because of a print-inside-list-comprehension and use of fstrings.
sed would work here, and you should learn about regular expressions (or regexs). If you aren't familiar with them. sed is a tool to edit files generally based on a pattern matching syntax, and regular expression are how those patterns are described. This might be a good start to learn regexs: https://regexone.com/ They'll likely become a very useful tool as you get more into bioinformatics.
The solution is to tell sed to "modify any line that doesn't start with a >, and any other line should be truncated to keep only the last 1000 characters. Since you know this is a FASTA file, you can assume than any line that starts with > should be left untouched. Then, your regex needs to do the truncation.
Here's an example to illustrate: I created a dummy file with two lines (one header row, and one row with 999 As and 1001 Bs), so the expected output is ABBBB...(such that there's 999 B's.) The sed command is:
$ sed 's|^[^>].*\(.\{1000\}\)|\1|g' file.txt
>fasta_header
ABBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
This may be really hard to understand and parse out on the first go and naturally regexes are not intuitive at first, you really just have to learn the commands and syntax. Don't get discouraged if you don't understand it all, just keep in mind that it's something that's possible to do with tools like sed and regexes.
Detail about the command, if you're interested: Sed will evaluate the expression line by line in the text file. The syntax for a replacement is `sed 's|pattern|replacement|g'. There's other types of sed command too.
The regex, |^[^>].*\(.\{1000\}\)|\1| means matchin from the begining of the file (^) any character that isn't > (that's [^>], it means "not >" essentially.), for an unlimited amount of characters (.*) up until any other character times 1000, .*\(.\{1000\}\). Because you used .*, sed will use greedy matching and go all the way until the end of the string, leaving you the last 1000 characters.
The \ here aren't actually part of the regex and are just there so sed doesn't interpret the ( ) { } as part of the expression, but rather as a regular expression group ( ) and a repeat character {NUMBER}. Really it should be (.{1000}) which means a group of any 1000 characters.
The group ( ) is a special regular expression construct that let's you re-use part's you've matched of your regex. Since it's the first (and only) group in the expression, the second part of the sed command (\1) means replace whatever the regular expression matched (if there's a match, so not applicable to headers) to whatever is inside the group (which should be the last 1000 characters).
Log in to answer this question.
See this post, which uses
awkto do some similar: https://stackoverflow.com/questions/30264809/trim-first-n-bases-in-multi-fasta-file-with-awk-and-print-with-max-width-formatThe command in the last code block should be somewhat easy to modify. Instead of using the left_trim argument to trim the first 1000 bases as in the example, you could modify the script to keep the last 1000 bases (something like replacing
start=left_trim+1;forstrat=length(rec) - left_trim + 1). Try it out and test on some examples.I could do it! I'm very thankful for your help.
Could you please paste here the explanation that you did about sed syntax?
I removed it because it didn't handle sequences that spanned multiple lines, which is common for FASTA files.
sedworks line by line, so awk was more appropriate in that use case. But the command wascat file.txt | sed 's|^[^>].*\(.{1000}\)|\1|g, if you're intrested. Briefly,sed 's|PATTERN|REPLACEMENT|g'will match lines that havePATTERNand substitute it withREPLACEMENT. The pattern is a regular expression^[^>].*\(.{1000}\), which means from the a line that begins (^) with anything except a>(^[^>]) followed by anything of unlimited length.*until you can't match anything but 1000 characters of anything (.{1000}, and see greedy matching which applies here.). The last part is wrapped in parenthesis (\(.{1000}\), or(GROUP)) whereGROUPis a part of the regex you want to keep for later. TheREPLACEMENTpart just says "replace with group 1 " (\1). The backslashes before{ } ( )are there to say you don't want to match those characters, instead they are part of the regular expression syntax. I'd highly recommend learning more aboutsedand regexes for bioinformatics, they're really useful and widely used.seqkit subseq supports this, if you want a fast solution.
If you want to learn programming, write some Python scripts using Biopython.