I've downloaded the RefSeq transcript sequences (FASTA format) for the GRC38 patch 14 build of the human genome and would now like to extract the gene symbols (only) in order from the headers. I'm confident I can do this using grep and regular expressions--the trick will be coming up with an appropriate pattern. (The gene symbol appears to be surrounded by the last pair of parentheses in a header, but I'm not sure this is absolutely true.) However, I'm wondering if there is a more elegant way to do this--e.g., via the NCBI EDirect utilities or SeqKit--that I'm overlooking. (I have a list of the accession.version numbers.) Any thoughts on this would be greatly appreciated.
$\begingroup$
$\endgroup$
3
-
5$\begingroup$ It would be helpful to see a sample of the file or a direct link to the file in question. Is there a reason that you want to do this from the FASTA header (which is probably not so strictly curated) rather than from a more traditional annotation file format such as GFF3? $\endgroup$– Maximilian PressCommented Jun 14, 2022 at 18:39
-
2$\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$– Community BotCommented Jun 14, 2022 at 19:27
-
$\begingroup$ We can't help you parse data if you don't show us what you need to parse. Please edit your question and show us a few lines (preferably just fasta header lines) of your file and the output you would want from that example. $\endgroup$– terdon ♦Commented Jun 21, 2022 at 10:27
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
2
We can't be sure without seeing your file, but if what you really want is the text between the last parentheses on lines that start with >, then you can do this with GNU grep (the default on Linux) with:
grep -oP '^>.*\(\K[^)]+' file.fa
The -o means "print only the matching portion of the line" and the -P activates Perl Compatible Regular Expressions (PCREs). The regex looks for a > at the beginning of the line (^>), then as many characters as possible until the last parenthesis (.*\(), and then forgets everything matched up to this point with the \K. This means that the -o will only see things after the \K, a PCRE feature. So finally, we look for the longest stretch of one or more non-) characters.
If you don't have GNU grep, you can do the same basic thing with perl or sed:
sed -n 's/^>.*(\([^)]*\).*/\1/p' file.fa
perl -lne '/^>.*\(([^)]+)/ && print $1' file.fa
$\begingroup$
$\endgroup$
4
There are probably better ways of doing this with Python or another tool, but here's a quick method that works in bash (tested on Linux) assuming that the Gene Symbols are contained in parentheses, i.e. "(Gene Symbol)":
pat='\(([^)]+)\)'
while read line; do
if [[ "$line" =~ $pat ]]; then
echo ${BASH_REMATCH[1]};
fi;
done < human.1.rna.fna
-
$\begingroup$ Thanks. Unfortunately, there can be more than one parenthesized expression in the header, so this pattern is too permissive. In every record I've looked at, the gene symbol is in the last parenthesized expression, but I don't know if that is guaranteed to be true. $\endgroup$ Commented Jun 20, 2022 at 14:53
-
$\begingroup$ Could you maybe use another file provided by the RefSeq ftp such as this one? It may have more genes than in your FASTA file, I'm assuming not fewer, but you could always do a quick check using the regex output as a comparison by sorting both files, and performing a
diffcommand in Linux. $\endgroup$– helrickCommented Jun 20, 2022 at 17:09 -
$\begingroup$ You really, really don't want to do this in the shell. It is going to be orders of magnitude slower than anything else and is very hard to get right. For example, your entire loop can be done faster and simpler with
grep '\([^)]*\) human.1.rna.fna. The loop only makes it slower and more fragile. $\endgroup$– terdon ♦Commented Jun 21, 2022 at 10:26 -
$\begingroup$ Relevant: Why is using a shell loop to process text considered bad practice? $\endgroup$– terdon ♦Commented Jun 21, 2022 at 10:37