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$
3
  • $\begingroup$ Could you provide the lines and the context around the lines? $\endgroup$
    – Bioathlete
    Feb 20, 2019 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$ Feb 20, 2019 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, 2019 at 9:12

4 Answers 4

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.

$\endgroup$
2
  • $\begingroup$ So, you don't want to get rid of them and process the file as a single object. $\endgroup$
    – M__
    Feb 21, 2019 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$ Feb 21, 2019 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")

$\endgroup$
2
  • $\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$
    – M__
    Feb 21, 2019 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$
    – M__
    Feb 22, 2019 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
$\endgroup$
5
  • $\begingroup$ Thanks. What does your awk command do? $\endgroup$ Feb 21, 2019 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, 2019 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, 2019 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$ Feb 21, 2019 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, 2019 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
$\endgroup$
9
  • $\begingroup$ what is your sed command doing? There is no output and no change to the file when I run it. $\endgroup$ Feb 20, 2019 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$
    – M__
    Feb 21, 2019 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, 2019 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$
    – M__
    Feb 21, 2019 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$
    – M__
    Feb 21, 2019 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.