3
$\begingroup$

I downloaded the annotation of the C. elegans genome in GFF3 format from Ensembl. I typed the following command, hoping to get the header of the file (lines starting with #).

grep '^#' Caenorhabditis_elegans.WBcel235.95.gff3

This returned 46798 lines, many many more than the dozen or so header lines. After searching the file for the # character, I noticed many lines that contain only three hashtag (#) symbols.

screenshot showing hashtags in GFF3 file

Why are these lines there?

$\endgroup$
  • $\begingroup$ Could you provide the lines and the context around the lines? $\endgroup$ – Bioathlete Feb 20 '19 at 23:00
  • $\begingroup$ I added a screenshot. There are too many lines to show them all here. There is no constant interval between the ### lines. $\endgroup$ – charlesdarwin Feb 20 '19 at 23:45
  • 3
    $\begingroup$ Next time, please don't post images of text. Just copy the text directly into your question and use the formatting tools. images are harder to search, make the page heavier to load (important for users with slow connection) and cannot be copied (I can't extract your file to demonstrate a command that solves your problem). $\endgroup$ – terdon Feb 21 '19 at 9:12
4
$\begingroup$

There are some good answers so far, but I don't think any of them fully communicate the significance of the ### directive. The GFF3 specification states:

This directive (three # signs in a row) indicates that all forward references to feature IDs that have been seen to this point have been resolved. After seeing this directive, a program that is processing the file serially can close off any open objects that it has created and return them, thereby allowing iterative access to the file. Otherwise, software cannot know that a feature has been fully populated by its subfeatures until the end of the file has been reached. It is recommended that complex features, such as the canonical gene, be terminated with the ### notation.

In other words, without the ### directive you have to load the entire GFF3 file into memory to make sure you've inferred all feature relationships completely (as encoded by ID and Parent attributes). With the ### directive, it's possible to process the annotations in smaller chunks, with a MUCH smaller peak memory usage.

| improve this answer | |
$\endgroup$
  • $\begingroup$ So, you don't want to get rid of them and process the file as a single object. $\endgroup$ – Michael Feb 21 '19 at 17:04
  • $\begingroup$ @MichaelG Right-o! Note, that having ### directives only makes it possible to process GFF3 files efficiently. It is not a guarantee that any particular tool will utilize these directives. $\endgroup$ – Daniel Standage Feb 21 '19 at 17:42
4
$\begingroup$

The three # are used for splitting group of features that belong together, e.g. a transcript and it's exons. Sometimes you see with a blank line instead of the #.

To grep just the header, which has 1 or 2 # from the line start on, you can use extended regular expression:

$ grep -E "^#{1,2}[^#]" Caenorhabditis_elegans.WBcel235.95.gff3

Which means: "grep all lines that started with 1 or 2 # and are not followed by #"

If you like to create a new gff3 to get rid of this separator lines you can do this:

$ grep -v "^###" Caenorhabditis_elegans.WBcel235.95.gff3 

The -v prints all line that doesn't match the pattern. (I guess this is what @Michael mean by "negator")

| improve this answer | |
$\endgroup$
  • $\begingroup$ I know "negator" is made up ... The technical term is "negation operator", (en.wikipedia.org/wiki/Negation). awk is more powerful than sed, as finswimmer demonstrates. ... original post edited $\endgroup$ – Michael Feb 21 '19 at 15:49
  • $\begingroup$ I discussed with @terdon. I didn't realise sed/grep had negation operator [^#], I thought this was regex such as Perl.. ... so .... grep "^#.[^#]" ./Caenorhabditis_elegans.WBcel235.95.gff3 ... is another solution $\endgroup$ – Michael Feb 22 '19 at 15:48
3
$\begingroup$

As others have stated, those are just there to separate the entries for easier parsing. They enable you to do nifty tricks like:

$ awk -v RF='###' '/Y74C9A.3/' Caenorhabditis_elegans.WBcel235.95.gff3 
I   WormBase    mRNA    4116    10230   .   -   .   ID=transcript:Y74C9A.3;Parent=gene:WBGene00022277;Name=Y74C9A.3;biotype=protein_coding;transcript_id=Y74C9A.3
I   WormBase    three_prime_UTR 4116    4220    .   -   .   Parent=transcript:Y74C9A.3
I   WormBase    exon    4116    4358    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e5;constitutive=1;ensembl_end_phase=-1;ensembl_phase=0;exon_id=Y74C9A.3.e5;rank=5;version=1
I   WormBase    CDS 4221    4358    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    5195    5296    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e4;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e4;rank=4;version=1
I   WormBase    CDS 5195    5296    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    6037    6327    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e3;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e3;rank=3;version=1
I   WormBase    CDS 6037    6327    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    9727    9846    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e2;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e2;rank=2;version=1
I   WormBase    CDS 9727    9846    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    CDS 10095   10148   .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    10095   10230   .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e1;constitutive=1;ensembl_end_phase=0;ensembl_phase=-1;exon_id=Y74C9A.3.e1;rank=1;version=1
I   WormBase    five_prime_UTR  10149   10230   .   -   .   Parent=transcript:Y74C9A.3

To extract the headers, you just need the lines that start with exactly 2 #:

$ grep '^##[^#]' Caenorhabditis_elegans.WBcel235.95.gff3 
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942

Or, to also include the comments, lines that start with a #, then either a # or a ! and then a non-#:

$ grep '^#[#!][^#]' Caenorhabditis_elegans.WBcel235.95.gff3 
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942
#!genome-build WormBase WBcel235
#!genome-version WBcel235
#!genome-date 2012-12
#!genome-build-accession NCBI:GCA_000002985.3
#!genebuild-last-updated 2017-06

Or just a #, any character and then a non-#:

$ grep '^#.[^#]' Caenorhabditis_elegans.WBcel235.95.gff3 
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942
#!genome-build WormBase WBcel235
#!genome-version WBcel235
#!genome-date 2012-12
#!genome-build-accession NCBI:GCA_000002985.3
#!genebuild-last-updated 2017-06
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks. What does your awk command do? $\endgroup$ – charlesdarwin Feb 21 '19 at 12:24
  • $\begingroup$ @charlesdarwin it prints the output I show in my answer :) Basically, it sets the string ### as the record ("line") separator and then print the record that matches Y74C9A.3. The default action when something evaluates to true in awk is to print the line. So, /Y74C9A.3/ will be true on lines that contain Y74C9A.3, and since we have used ### as the line separator (record separator) it will print the entire entry for Y74C9A.3. $\endgroup$ – terdon Feb 21 '19 at 12:37
  • $\begingroup$ @charlesdarwin why did you add the space after the -v? Do you know of an awk implementation that requires it? I fear that makes it look like a "verbose" flag. $\endgroup$ – terdon Feb 21 '19 at 17:46
  • $\begingroup$ my awk complained without the space. awk: invalid -v option. I am on macOS X high sierra. using awk version 20070501. crazy number there! is it a date? $\endgroup$ – charlesdarwin Feb 21 '19 at 21:24
  • 1
    $\begingroup$ @charlesdarwin ah, yes, I tested on GNU awk. Thanks! And yes, that looks like the 1st of May, 2007. $\endgroup$ – terdon Feb 21 '19 at 21:27
-1
$\begingroup$

Hmmm, I dunno. If I was doing this I'd want an easy way to parse the file and ### would be a good marker. Anyway just pop a space character after the #. Using sed ...

sed -n '/^# /p' < ./Caenorhabditis_elegans.WBcel235.95.gff3

That should remove the ### lines

Also head -20 and sed -n '/1,20/p' ./Caenorhabditis_elegans.WBcel235.95.gff3 would get around the problem with a bit of tweaking.


The grep verion is (removing unnecessary jargon)

grep '^# ' Caenorhabditis_elegans.WBcel235.95.gff3

This is the same command,

grep '^#\s' Caenorhabditis_elegans.WBcel235.95.gff3

The point behind the code is there is no negation operator, in sed (I assume grep is the same). In Perl you would simple [^#]{3} (need checking) in your regex (remove all ###). In sed (or grep) you can exclude them using the ^ (beginning of a line) and suffixed with a "space" with only a single character in the middle all 3 letter patterns are removed.


Another approach without the space character is

sed -n '/^###/d; /^#/p' ./Caenorhabditis_elegans.WBcel235.95.gff3

Albeit this may require,

sed -i '/^###/d; /^#/p' ./file
| improve this answer | |
$\endgroup$
  • $\begingroup$ what is your sed command doing? There is no output and no change to the file when I run it. $\endgroup$ – charlesdarwin Feb 20 '19 at 23:48
  • $\begingroup$ I've edited the post to include the command for your grep approach and included an explanation. It is worth learning sed BTW, it has some nice tricks, but it is not as important as awk. $\endgroup$ – Michael Feb 21 '19 at 0:02
  • $\begingroup$ Your sed command does nothing. I think you meant to use -i there, but I don't really see how. There are no lines matching ^# in the file. To remove the ### lines, you would want something like sed -i '/^###/d' file.gff. And you can negate character classes in both sed and grep just as you do in perl: [^#] means "not #". This is a general regex feature available in all common regex flavors. It isn't a PCRE thing (although GNU grep can deal with PCRE as well, just use the -P flag). But note that [^###] does not mean "not ###", it just means "not #". $\endgroup$ – terdon Feb 21 '19 at 9:31
  • $\begingroup$ sed -i '/^###/d' would not work BTW, becaue sed -n '/^###/d' does not. sed '/^###/d' is an improvement. One way to get it working would be >sed -n '/^##/d; /^#/p' ./file ... I've edited the post accordingly $\endgroup$ – Michael Feb 21 '19 at 15:50
  • $\begingroup$ I suspect sed -n has been deprecated with implementations of sed -i ... sed -n is discussed here stackoverflow.com/questions/6255796/… … ...... It appears sed -i is not supported in all sed implementations, stackoverflow.com/questions/18527365/what-does-sed-i-option-do .... SED (-n): is direct from O'Reilly's famous "sed & awk" Doughty & Robbins, e.g p71 ... sed -n works on OS X and Susie [bash and tcsh]. PERL: [^###] is identical to [^#] but [^#]{3} I think would resolve the issue. $\endgroup$ – Michael Feb 21 '19 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.