text:
How to indicate the END of a haplotype block in VCF?
In VCF I know how to indicate that two genotypes are in the same phase by using consecutive "0|1" and "1|0" genotype fields, for example. However, how do I deal with the case that the first two genotypes are in phase and the second 2 are in phase, but there is no phase between the 2nd and 3rd?
According to the format I think I need to write:
chr1 100 . A T . . . GT 1|0
chr1 200 . A T . . . GT 0|1
chr1 300 . A T . . . GT 1|0
chr1 400 . A T . . . GT 0|1

but that implies that the SNP at 200 and the SNP at 300 are in phase, when in fact they are not.
How to indicate the END of a haplotype block in VCF?

According to the VCF specification, haplotype blocks are defined by the "Phase Set" [PS] tag:

PS : phase set. A phase set is defined as a set of phased genotypes to which this genotype belongs. Phased genotypes for an individual that are on the same chromosome and have the same PS value are in the same phased set. A phase set specifies multi-marker haplotypes for the phased genotypes in the set. All phased genotypes that do not contain a PS subfield are assumed to belong to the same phased set. If the genotype in the GT field is unphased, the corresponding PS field is ignored. The recommended convention is to use the position of the first variant in the set as the PS identifier (although this is not required). (Non-negative 32-bit Integer)

A good example of this can be found on a 10X Genomics support page about phasing:
CHROM   Pos REF ALT GT  PS
chr1    1000    A   C   0|1 1000
chr1    1010    T   G   1|0 1000
chr1    2000    C   T   0|1 2000
chr1    2005    T   G   0/1 2000
chr1    2008    G   C   0|1 2000

In this example we have two phase blocks, denoted by PS=1000 and PS=2000. PS=1000 spans position 1000-1010, and PS=2000 spans position 2000-2008. In PS=1000, haplotype 1 contains the REF A allele at position 1000, and the ALT G allele as position 1010, while haplotype 2 contains the ALT C allele at position 1000 and the REF T allele at position 1010.

Note that this phase set identifier needs to be a single numeric value, which means that sub-blocks and overlapping blocks cannot be defined by this field.

Two options that immediately come to mind, depending on what you hope to accomplish down the pipeline: you might simply mark the unphased sites with a "/" instead of "|" to indicate they are not phased; or, if you have multiple haplotype blocks that are all independent of each other, you could use the chromosome field to identify each unique haplotype (e.g. chr1h1, chr1h2, chr2h1, chr2h2, etc).

There is a new idea of reporting the "phased haplotype/genotype block" with unique block ID.
phaser uses "PI" (unique phase-index) and "PG" (phased-genotype).
GATK is also using this idea of "unique-block" and "phased-genotype"; See the "HP" tag.
Personally, I feel the representation of "phaser" is way better and comprehensive.
I would represent your data as:
chr1 100 . A T . . . GT:PI:PG    1/0:3:1|0
chr1 200 . A T . . . GT:PI:PG    0/1:3:0|1
chr1 300 . A T . . . GT:PI:PG    1/0:4:1|0
chr1 400 . A T . . . GT:PI:PG    0/1:4:1|0

where first two lines have same "PI" (i.e 3) stating they are phased with each other. Same with the later two lines which have another "PI" (i.e 4).


text:
What is a good pipeline for using public domain exomes as controls?
I'm currently attempting association analysis with an extremely small set of patient exomes (n=10), with no control or parental exomes available. Downloading the ExAC VCF of variant sites (http://exac.broadinstitute.org/downloads) or the 1000G integrated call sets (http://ftp.1000genomes.ebi.ac.uk/) and combining this with our pooled patient VCFs has not been successful (I suspect the approach of attempting to merge such large VCFs generated from different pipelines is rather naive).
Looking at the primary literature, I have gathered it should be possible to use these resources to help increase statistical power for our analysis. My question is how do I take these large .vcfs with many samples and successfully merge them to our patient .vcfs, such that the combined VCF can be used downstream to run analysis packages? (PODKAT, PLINK, etc.)

I gather that you want to use the background SNP frequency as a prior for input to your SNP calling algorithm?
I'm not sure of a canned algorithm for doing this, but a quick google shows up some promising links:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3593722/
https://software.broadinstitute.org/gatk/documentation/article.php?id=4723
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0489-0

Thanks for the response - to clarify, I need additional details on how to correctly assign my priors using these public resources to jointly analyse multiple samples (our patient exomes and the population controls). The first link touches on this, but I'm unsure how to actually implement this approach.
The third link would be useful, except our coverage is fine and instead we lack the sample size they have by several orders of magnitude (gulp)!

Yeah, I don't know of any tools that do this using existing data. What you /can/ do is re-analyze your data in the context of the other data in a single multi-sample call. This seems like a hack TBH, but should work.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
RATT works on example bacterial sequence and other bacterial genome but not on C. elegans genome and annotation
I am trying to use a tool called RATT described in this paper and available here. I ran the example using the following command:
$RATT_HOME/start.ratt.sh embl F11.fasta F11 Strain
embl is a folder containing the file Tb_H37Rv.embl, which is an annotation file in EMBL format of the Mycobacterium tuberculosis genome. F11.fasta contains the genome in FASTA sequence of the strain F11 of Mycobacterium tuberculosis. This example provided by the RATT developers runs smoothly:
Of the reference chromosome Tb_H37Rv    2.66 per cent   has no synteny with the query
Of the query chromosome gi_148821191_ref_NC_009565.1_    3.03 per cent  has no synteny with the reference

...
Please see the file F11.gi_148821191_ref_NC_009565.1_.Report.txt for reporting the errors and the file F11.gi_148821191_ref_NC_009565.1_.Report.gff for reporting the error in Artemis or other sequence viewer.
************************
If you want to start artemis on this replicon:
art F11.gi_148821191_ref_NC_009565.1_.final.embl + F11.gi_148821191_ref_NC_009565.1_.Report.gff  + Query/F11.gi_148821191_ref_NC_009565.1_.Mutations.gff

Then, I tried it on the genome of E. coli and it worked:
$RATT_HOME/start.ratt.sh embl bact1.fa bact1 Species

I used the embl file from Mycobacterium so not much annotation was transferred due to the large genetic differences between Mycobacterium and E. coli.
Then, I tried to transfer the annotation from the reference genome of C. elegans to a new assembly I have generated. The way I generated the embl annotation file is by running seqret on the fasta file of C. elegans and its corresponding annotation file found in GTF format:
seqret -sequence c_elegans.PRJNA13758.WS263.genomic.fa -feature -fformat gff -fopenfile c_elegans.PRJNA13758.WS263.annotations.gff3 -osformat embl -auto -outseq c_elegans.PRJNA13758.WS263.embl

I got the following error:
RATT_HOME/start.ratt.sh embl scaffolds.fa scaffolds Strain

ERROR: Reference input does not match delta file
Assertion failed: (Len - Coord + 1 > 0), function revC, file show-coords.cc, line 1118.
/Users/user/Documents/project/output/software/ratt-code/start.ratt.sh: line 73: 85325 Abort trap: 6           show-coords -clHT $name.delta > $name.coords
Assertion failed: (Len - Coord + 1 > 0), function revC, file show-coords.cc, line 1118.
/Users/user/Documents/project/output/software/ratt-code/start.ratt.sh: line 73: 85326 Abort trap: 6           show-coords -clHT $name.filter.delta > $name.filter.coords
Use of uninitialized value $stats in print at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 82.
Nucmer is done. Now transfer the annotation.
perl /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl embl nucmer.scaffolds.snp nucmer.scaffolds.filter.coords scaffolds
working on c_elegans.PRJNA13758.WS263
Use of uninitialized value in split at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 1140.
Can't use an undefined value as an ARRAY reference at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 229.
6 Sequences where generated
Done.
Nucmer is done. Now Correct the annotation for chromosome .
work on I
If you want to start artemis on this replicon:
art scaffolds.I.final.embl + scaffolds.I.Report.gff  + Query/scaffolds.I.Mutations.gff
work on II
If you want to start artemis on this replicon:
art scaffolds.II.final.embl + scaffolds.II.Report.gff  + Query/scaffolds.II.Mutations.gff
work on III
If you want to start artemis on this replicon:
art scaffolds.III.final.embl + scaffolds.III.Report.gff  + Query/scaffolds.III.Mutations.gff
work on IV
If you want to start artemis on this replicon:
art scaffolds.IV.final.embl + scaffolds.IV.Report.gff  + Query/scaffolds.IV.Mutations.gff
work on V
If you want to start artemis on this replicon:
art scaffolds.V.final.embl + scaffolds.V.Report.gff  + Query/scaffolds.V.Mutations.gff
work on X
If you want to start artemis on this replicon:
art scaffolds.X.final.embl + scaffolds.X.Report.gff  + Query/scaffolds.X.Mutations.gff
rm: scaffolds*embl.tmp.BBA.embl: No such file or directory

I downloaded the FASTA file here and the GTF file here.
I am working on Mac OSX.
Update 6th August 2018
I have also tried an alternative way of getting a C. elegans annotation in EMBL format. I downloaded it from the Ensembl FTP website:
wget -m ftp://ftp.ensembl.org/pub/release-91/embl/caenorhabditis_elegans/ --directory-prefix embl/
Then, I cat-ed the .dat files (one for each chromosome) into one ensembl.embl file and used that for the RATT command:
$RATT_HOME/start.ratt.sh embl scaffolds.fa scaffolds Strain

ERROR: Could not parse delta file, nucmer.scaffolds.delta
error no: 402
ERROR: Could not parse delta file, nucmer.scaffolds.filter.delta
error no: 402
ERROR: Could not parse delta file, nucmer.scaffolds.delta
error no: 402
ERROR: Could not parse delta file, nucmer.scaffolds.filter.delta
error no: 402
Use of uninitialized value $stats in print at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 82.
Nucmer is done. Now transfer the annotation.
perl /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl embl nucmer.scaffolds.snp nucmer.scaffolds.filter.coords scaffolds
working on ensembl
Use of uninitialized value in split at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 1140.
Can't use an undefined value as an ARRAY reference at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 229.
6 Sequences where generated
Done.
Nucmer is done. Now Correct the annotation for chromosome .
work on I
If you want to start artemis on this replicon:
art scaffolds.I.final.embl + scaffolds.I.Report.gff  + Query/scaffolds.I.Mutations.gff

...
work on X
If you want to start artemis on this replicon:
art scaffolds.X.final.embl + scaffolds.X.Report.gff  + Query/scaffolds.X.Mutations.gff
rm: scaffolds*embl.tmp.BBA.embl: No such file or directory

Indeed, the nucmer delta file is empty. If that is the only problem, maybe I can run nucmer separately but then how can I run RATT without running nucmer?

The developer of RATT told me to not concatenate the EMBL files because there is no need to do so. So, here is the output of RATT on the .embl files, one for each chromosome.
$RATT_HOME/start.ratt.sh embl scaffolds.fa scaffolds Strain

Of the reference chromosome Caenorhabditis_elegans.WBcel235.91.chromosome.I 13.76 per cent  has no synteny with the query

...
Of the query chromosome X    15.29 per cent has no synteny with the reference
Nucmer is done. Now transfer the annotation.
perl /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl embl nucmer.scaffolds.snp nucmer.scaffolds.filter.coords scaffolds
working on Caenorhabditis_elegans.WBcel235.91.chromosome.I
Argument "11123466)" isn't numeric in sort at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 985, <F> line 2547.
Argument "11123466)" isn't numeric in sort at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 985, <F> line 2555.
Argument "5568506)" isn't numeric in sort at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 985, <F> line 2627.
1474993 1480045 2526 1477519
14358883 14373963 7540 14366423
3204189 3216229 6020 3210209
3809476 3812106 1315 3810791
10232464 10233260 398 10232862
10491944 10506864 7460 10499404
11020347 11023081 1367 11021714
318712 361504 21396 340108
working on Caenorhabditis_elegans.WBcel235.91.chromosome.MtDNA
Can't use an undefined value as an ARRAY reference at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 229.
6 Sequences where generated
Done.
Nucmer is done. Now Correct the annotation for chromosome .
work on I
If you want to start artemis on this replicon:
art scaffolds.I.final.embl + scaffolds.I.Report.gff  + Query/scaffolds.I.Mutations.gff

...
work on X
If you want to start artemis on this replicon:
art scaffolds.X.final.embl + scaffolds.X.Report.gff  + Query/scaffolds.X.Mutations.gff
rm: scaffolds*embl.tmp.BBA.embl: No such file or directory

It is better than before because I get a reasonable delta file from nucmer. But the error Argument "11123466)" isn't numeric is a problem. It seems that is should be easy enough to solve this. It looks to me like the coordinates from the EMBL file are not parsed properly and the code stumbles on the coordiate followed by a ).

The seqret didn't give any error?  It seems like the expectation of ratt and the output of seqret differ. Maybe it needs a specific version or the problem is only on the RATT program

seqret did not give an error.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Can I export the SNPs between whole genome alignments using the command line?
I am using Mauve to align two whole genomes. I use the following command to get the alignment in xmfa format:
/Applications/Mauve.app/Contents/MacOS/progressiveMauve ref1.fa ref2.fa --output alignment.xmfa
I then open the xmfa file in the Mauve GUI and export the SNPs. Is there a way of doing this using a command? I have tried to follow the procedure described here but the link to snapshot releases of Mauve does not work anymore.
I could not find anything in the Mauve command line help:
/Applications/Mauve.app/Contents/MacOS/progressiveMauve --help
nor on the Mauve website.
I tried the procedure described by Konrad Rudolph. I wrote:
MAUVE_DIR=/Applications/Mauve.app/Contents/Java
in the terminal. Is this correct? Don't I need an export and quotation marks to create a bash variable? I tried adding an export and using quotation marks for the file but the error below remained the same.
followed by:
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
followed by:
java org.gel.mauve.analysis.SnpExporter -f alignment.xmfa -o aln.snps
and it seems to work.

It works for me with the download from download page. The download archive contains Mauve.jar so the command line in the answer you link to should work.
However, you might need to include all the jar files (not just Mauve.jar) in the classpath for this to command work, though. For instance, the following works for me:
MAUVE_DIR=‹the-dir-containing-Mauve.jar›
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
java org.gel.mauve.analysis.SnpExporter -f aln.xmfa -o aln.snps

… there’s probably a smarter way to make this work but I don’t use Java enough to know how.

Should I download this file: http://darlinglab.org/mauve/snapshots/2015/2015-02-25/MacOS/Mauve-snapshot_2015-02-25.dmg ? But this is to install mauve, right? Where can I find the Mauve.jar file?

@charlesdarwin macOS applications are regular directories. Open the DMG and go inside the Mauve.app dir. Inside that directory you’ll find /Contents/Java/Mauve.jar; to display the contents in Finder, use the instructions here: https://apple.stackexchange.com/a/95942/2392 — but in the terminal you can simply cd into the directory.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
What is the `table` parameter in the UCSU ldHgGene tool?
I want to lift over the annotation from E. coli release 29 to E. coli release 42. I have followed the steps for an annotation liftover mentioned here. I got the chain file for the final step and ran:
liftOver -gff Escherichia_coli.GCA_000967155.1.gtf Chromosome.chain conversions.gtf unmapped.gtf
However, I get the following warning:
WARNING: -gff is not recommended.
Use 'ldHgGene -out=<file.gp>' and then 'liftOver -genePred <file.gp>'

I would like to heed the warning. However, I don't know how to use the ldHgGene tool. Here is its help page:
ldHgGene - load database with gene predictions from a gff file.
usage:
     ldHgGene database table file(s).gff
options:
     -bin         Add bin column (now the default)
     -nobin       don't add binning (you probably don't want this)
     -exon=type   Sets type field for exons to specific value
     -oldTable    Don't overwrite what's already in table
     -noncoding   Forces whole prediction to be UTR
     -gtf         input is GTF, stop codon is not in CDS
     -predTab     input is already in genePredTab format
     -requireCDS  discard genes that don't have CDS annotation
     -out=gpfile  write output, in genePred format, instead of loading
                  table. Database is ignored.
     -genePredExt create a extended genePred, including frame
                  information and gene name
     -impliedStopAfterCds - implied stop codon in GFF/GTF after CDS

From what I understand, I am trying to give the liftOver tool the annotation in GenePred format instead of GTF format. So I need to convert the GTF format into a GenePred format and that's what the ldHgGene tool is for. So the command should be something like:
ldHgGene Escherichia_coli.GCA_000967155.1.gtf -out=Escherichia_coli.GCA_000967155.1.gp

but this returns nothing. The -out option does mean that the database parameter is not needed but the tool still seems to need the parameter table.
What is the table parameter?

I stumbled across this issue after getting the same warning message from liftOver. By using the -out option of lgHgGene, the database parameter is ignored but must still be provided. The table parameter should be the second column from your GTF file. Note that you must also specify the -gtf option (otherwise it expects GFF format). In my case, the GTF file looks like:
chr1    mm10_knownGene  exon    3073253 3074322 0.000000        +       .       gene_id "ENSMUST00000193812.1"; transcript_id "ENSMUST00000193812.1";
chr1    mm10_knownGene  exon    3102016 3102125 0.000000        +       .       gene_id "ENSMUST00000082908.1"; transcript_id "ENSMUST00000082908.1";
chr1    mm10_knownGene  exon    3206523 3207317 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3213439 3215632 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3205901 3207317 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
chr1    mm10_knownGene  exon    3213609 3216344 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
...

By default, ldHgGene outputs a leading 'bin' column which won't work with liftOver; this can be avoided with the -nobin option. So my command was:
ldHgGene -gtf -nobin -out=mm10.GENCODE_VM20.gp ignored mm10_knownGene mm10.GENCODE_VM20.gtf

The resulting genePred output can be fed directly into liftOver -genePred:
ENSMUST00000193812.1    chr1    +       3073252 3074322 3074322 3074322 1       3073252,        3074322,
ENSMUST00000082908.1    chr1    +       3102015 3102125 3102125 3102125 1       3102015,        3102125,
ENSMUST00000162897.1    chr1    -       3205900 3216344 3216344 3216344 2       3205900,3213608,        3207317,3216344,
ENSMUST00000159265.1    chr1    -       3206522 3215632 3215632 3215632 2       3206522,3213438,        3207317,3215632,

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to use busco installed through conda?
I have installed busco on my Mac computer using conda install --channel bioconda busco=3.0.2 into an environment called busco3_0_2. However, when I type busco, the command is not found.
Anyone knows what the name of the binary is? Usually when I install something with conda, the name of the command is the same as the name of the package. In this case, it does not seem to be so.
Note that the package installation finished successfully.

From the README, it looks like the command-line entry point is not busco but run_BUSCO.py. This seems to agree with the bioconda recipe, where run_BUSCO.py and run_busco are listed among the commands.

Cheers. I found run_BUSCO.py and run_busco in ~/anaconda/envs/busco3_0_2/bin.

The new version of busco (v4) has the user entry point called busco (it's still a python script). BUSCO v4 is available at conda, note it's a very fresh update.


text:
Purpose of ### (3 consecutive pound signs / hashtags / octothorps) in GFF3
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.

Why are these lines there?

Could you provide the lines and the context around the lines?

I added a screenshot. There are too many lines to show them all here. There is no constant interval between the ### lines.

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).

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.

So, you don't want to get rid of them and process the file as a single object.

@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.

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")

I know "negator" is made up ... The technical term is "negation operator", (https://en.wikipedia.org/wiki/Negation). awk is more powerful than sed, as finswimmer demonstrates. ... original post edited

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

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

@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.

@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.

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?

@charlesdarwin ah, yes, I tested on GNU awk. Thanks! And yes, that looks like the 1st of May, 2007.

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

what is your sed command doing? There is no output and no change to the file when I run it.

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.

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 #".

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

I suspect sed -n has been deprecated with implementations of sed -i ... sed -n is discussed here https://stackoverflow.com/questions/6255796/how-the-n-command-works-in-sed … ...... It appears sed -i is not supported in all sed implementations, https://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.

No, sed -n is in no way deprecated, it's defined by POSIX.  And it has nothing to do with -i ("make the changes in the original file"), -n means "don't print unless I tell you to". And -i it is supported (with minor format differences) in all sed implementations you are likely to use as a bionformatician (GNU sed on Linux, BSD sed on BSD and OSX, even busybox sed on embedded systems). Anyway, sed -i '/^###/d'  works just fine, try it out! It's just that it edits the original file, so you don't see any output.

By the way,  sed -n '/^###/d' also "works", it just doesn't do anything useful. That means "delete lines starting with ### but in any case, don't ever print anything". Perhaps you meant something like sed -n '/^###/!p'.

There must by different versions of sed. I do not have any negation operators on mine - this was always my understanding of its weakness, ie your '/^###/!p'. I do agree your /^###/d' can be stacked into a /^#/p command and act as a work-around for negation, which I'd never thought about.

Let us continue this discussion in chat.


text:
bcftools consensus not working when using -H flag
Sorry, I have been trying to work this out desparately all day. Biostars couldn't find any answer so I am cross posting here.
I am trying to get both the left and right reference genome using a sample name in phase VCF file.
This works:
bcftools consensus -f lyrata_chr01-short.fasta phasedVCF-short.vcf.gz -s ms01e -H 1 > ms01e-left.fa

When I am trying to make the new Ref Sequence by imputing genotype from the right haplotype, this doesn't work.
bcftools consensus -f lyrata_chr01-short.fasta phasedVCF-short.vcf.gz -s ms01e -H 2 > ms01e-right.fa
Broken VCF, too few alts at 1:141

What is the issue ?
Here is my VCF data:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  MA605   MA611   MA622   MA625   MA629   Ncm8    Sp154   Sp164   Sp21    Sp3 Sp76    SpNor33 ms01e   ms02g   ms03g   ms04h
1   141 .   C   T   .   .   .   GT:PG:PG_al:PI  .   .   .   0|1:0|1:C|T:10763   .   .   .   .   .   .   .   .   .   .   .   .
1   3697    .   T   C   .   .   .   GT:PG:PG_al:PI  .   0|0:0|0:T|T:11398   .   0|1:0|1:T|C:8039    0|1:0|1:T|C:8202    0|1:0|1:T|C:4616    0|1:0|1:T|C:3938    .   0|1:0|1:T|C:5102    .   .   0|1:0|1:T|C:17340   1|0:1|0:C|T:1   .   .   .

The GT field is empty at POS 141 for SAMPLE ms01e, but in that situation it should by default pick a REF allele. Isn't it?
And, it only is the problem when doing -H 2 not -H 1.
Update: I have been trying to fix this issue by playing with -i and -e using the BCFtools manual. But, still cannot find the fix.

I don't know much about VCF but the first row seems to have less fields than the second. Is this normal/accepted format ? Also did you try the latest version (as per ) ?


text:
HISAT2 - are all multi mapped reads reported by default?
While going through the HISAT2 manual, one can see this:

Two alignments for the same individual read are "distinct" if they map
  the same read to different places. Specifically, we say that two
  alignments are distinct if there are no alignment positions where a
  particular read offset is aligned opposite a particular reference
  offset in both alignments with the same orientation. E.g. if the first
  alignment is in the forward orientation and aligns the read character
  at read offset 10 to the reference character at chromosome 3, offset
  3,445,245, and the second alignment is also in the forward orientation
  and also aligns the read character at read offset 10 to the reference
  character at chromosome 3, offset 3,445,245, they are not distinct
  alignments.
Two alignments for the same pair are distinct if either the mate 1s in
  the two paired-end alignments are distinct or the mate 2s in the two
  alignments are distinct or both.
Default mode: search for one or more alignments, report each

and then there is -k parameter (analogous to bowtie2 -K parameter)
which defines:

HISAT2 searches for up to N distinct, primary alignments for each
  read, where N equals the integer specified with the -k parameter.
  Primary alignments mean alignments whose alignment score is equal or
  higher than any other alignments. It is possible that multiple
  distinct alignments have the same score. That is, if -k 2 is
  specified, HISAT2 will search for at most 2 distinct alignments. The
  alignment score for a paired-end alignment equals the sum of the
  alignment scores of the individual mates. Each reported read or pair
  alignment beyond the first has the SAM 'secondary' bit (which equals
  256) set in its FLAGS field. See the SAM specification for details.

This is a bit confusing, as the first quote means that HISAT2 finds all distinct valid alignment and reports them by default, and the second quote says that when -k option is used, it will only report a subset of alignments which it would have reported by default. Am I correct in this inference or does one needs to use -k option to enable "multi-mapping"?

You said:

...the first quote means that HISAT2 finds all distinct valid alignment and reports them by default.

Not exactly. The first section you quoted doesn't say it will find and report all distinct alignments, but rather that it would find and report one or more distinct alignments. It's a subtle difference, but an important one.
That section describes how distinct mappings are defined, not how many are reported. The second section you quoted, as well as the documentation for the -k parameter, are all consistent with this interpretation.

-k <int>        It searches for at most <int> distinct, primary alignments for
                each read. Primary alignments mean alignments whose alignment
                score is equal or higher than any other alignments. The search
                terminates when it can't find more distinct valid alignments, or
                when it finds <int>, whichever happens first. The alignment
                score for a paired-end alignment equals the sum of the alignment
                scores of the individual mates. Each reported read or pair
                alignment beyond the first has the SAM 'secondary' bit (which
                equals 256) set in its FLAGS field. For reads that have more
                than <int> distinct, valid alignments, hisat2 does not guarantee
                that the <int> alignments reported are the best possible in
                terms of alignment score. Default: 5 (HFM) or 10 (HGFM)

                Note: HISAT2 is not designed with large values for -k in mind,
                and when aligning reads to long, repetitive genomes large -k can
                be very, very slow.

By default, 5 or 10 distinct mappings are reported (depending on the mode), but you can control this with the -k parameter.

Aha, thanks for clearing it up! Although I'm still a bit confused to why HISAT2 still has --secondary if setting up -k instructs it to search for multiple alignments n alignments regardless.

The -k parameter corresponds only to primary alignments: cases where distinct alignments have an identical high score. If an alignment has a less-than-optimal score, it is a secondary alignment.


text:
GATK CombineVariants complains the contig order in the VCF files
I have called variants on two strains of C. elegans separately. I now want to merge the VCF files into one using the following code:

Create a sequence dictionary of the reference sequence
Sort the VCF files with Picard
Merge the sorted VCF files using GATK

picard CreateSequenceDictionary \
    REFERENCE=c_elegans.PRJNA13758.WS263.genomic.fa \
    OUTPUT=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain1.vcf \
    OUTPUT=strain1sorted.vcf \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain2.vcf \
    OUTPUT=strain2sorted.vcf  \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

GATK --analysis_type CombineVariants \
    -R c_elegans.PRJNA13758.WS263.genomic.fa \
    --variant strain1sorted.vcf \
    --variant strain2sorted.vcf \
    -o all.vcf \
    -genotypeMergeOptions UNIQUIFY

The last command gives me the following error message:
ERROR MESSAGE: Input files variant and reference have incompatible contigs. Please see https://www.broadinstitute.org/gatk/guide/article?id=63 for more information. Error details: The contig order in variant and reference is not the same; to fix this please see: (https://www.broadinstitute.org/gatk/guide/article?id=1328),  which describes reordering contigs in BAM and VCF files..
##### ERROR   variant contigs = [I, II, III, IV, MtDNA, V, X]
##### ERROR   reference contigs = [I, II, III, IV, V, X, MtDNA]

But I have sorted the VCF files using Picard, so I don't know what else to do.
Your help is appreciated.
Update: I asked this question on the GATK forum on 4th May 2018 (today).

GATK cannot combine VCF files generated by Sniffles, the variant caller that I used to call structural variants (answer from the GATK team).

Answered in the GATK forum

Try regenerating the index files for your VCFs. Picard SortVcf doesn't do it for you iirc, so when GATK looks at the index files (before opening the VCFs themselves) it sees the old order and complains.

Hi, thanks for the answer, it's already helpful. However, It would be even nicer for community to write the actual answer and to provide the link just as the source (you can [edit] your answer).


text:
Batch alignment of inconsistently named Fastq files
I have numerous gzip-compressed paired-end Fastq files with ChIP-seq data that I would like to align with bowtie2. I confirmed that bowtie2 takes .gz files, although I do see this warning when I run for two samples:

Warning: gzbuffer added in zlib v1.2.3.5. Unable to change buffer size
  from default of 8192.

My files are named like so.
10_S10_L002_R1_001.fastq.gz  12_S12_L002_R1_001.fastq.gz  2_S2_L001_R1_001.fastq.gz  4_S4_L001_R1_001.fastq-003.gz  6_S6_L001_R1_001.fastq-002.gz  8_S8_L002_R1_001.fastq.gz
10_S10_L002_R2_001.fastq.gz  12_S12_L002_R2_001.fastq.gz  2_S2_L001_R2_001.fastq.gz  4_S4_L001_R2_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz  8_S8_L002_R2_001.fastq.gz
11_S11_L002_R1_001.fastq.gz  1_S1_L001_R1_001.fastq.gz    3_S3_L001_R1_001.fastq.gz  5_S5_L001_R1_001.fastq-002.gz  7_S7_L002_R1_001.fastq.gz      9_S9_L002_R1_001.fastq.gz
11_S11_L002_R2_001.fastq.gz  1_S1_L001_R2_001.fastq.gz    3_S3_L001_R2_001.fastq.gz  5_S5_L001_R2_001.fastq-003.gz  7_S7_L002_R2_001.fastq.gz      9_S9_L002_R2_001.fastq.gz

The bash commands I tried are below. When I run this test command:
for i in $(ls *.fastq*.gz | grep -v 6 | rev | cut -c 16- | rev | uniq); do
    echo $i;
done

...I get this output.
12_S12_L002_R1_001.fastq.gz  12_S12_L002_R2_001.fastq.gz

When I run this:
for i in $(ls *.fastq*.gz | rev | cut -c 16- | rev | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}R1_001.fastq.gz -2 ${i}R2_001.fastq.gz | samtools view -bS - >${i%}.bam
done

I do see one BAM for my paired end files. But as I see the second set of files which are named:
6_S6_L001_R1_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz

if I run the same script over them, I would get an error as the loop can't read the paired end file.
So the tedious way is I can make a separate folder for these sets of data with fastq-002.gz and fastq-003.gz and can run the script by tweaking the cut parameter.
But I would like to know what is the programmatic way of resolving the issue.

I strongly encourage you to use something like snakeMake rather than making a shell script. It's both much easier to ensure the proper files are paired and it scales to use in clusters and the cloud.

Can't you just rename the files?

@terdon that's how our collaborator gave us and i dont want to mess up with the file name yes i can if i dont find any other way..

" strongly encourage you to use something like snakeMake rather than making a shell script. It's both much easier to ensure the proper files are paired and it scales to use in clusters and the cloud." i read a bit about it never tried it , i need to learn and give it a try

It's good practice to treat files from collaborators as read-only, so not changing the filenames is a good idea. But if the files are inconsistently named, it's a very good idea to create consistently named symlinks to those files. Well organized filenames and directory structures will make it SO NICE AND EASY to automate batch data processing.

yes i did that finally

A very simple approach, but one which will work with the names you show, is to remove any occurrence of 00N (where N is any number) from right before the .gz. To do so, edit your script to:
for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do
  bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \
    -1 ${i}R1_001.fastq.gz \
    -2 ${i}R2_001.fastq.gz |
        samtools view -bS - >${i%}.bam"
done


text:
Autodetect max number of cores and pass as an argument in Nextflow
I am creating a pipeline in Nextflow. One step is creating a pangenome with Roary. Roary takes threads as an argument and if no number of threads is supplied as an argument it defaults to one.
Is there a way in Nextflow to pass the maximum number of available threads as an argument? At the moment the best I can come up with is defining an arbitrarily large number of threads as an argument which will cap out and run with the max number of threads. I would like to know, is there a 'cleaner' or more canonical way of doing this?
At the moment my process looks like this:

process roary {

    publishDir "${params.outdir}/roary", mode: 'copy'

    cpus 200

    input:
    file gff from gff.collect()

    output:
    file("*") into roary
    file("pan_genome_reference.fa") into pan_genome
    file("pan_genome_sequences/*") into alignment_files
    file("gene_presence_absence.Rtab") into gene_presence
    set file("*accessory*"),
        file("*.Rtab"),
        file("_*"),
        file("*.txt"),
        file("*csv") into roary_out

    script:
    """
    roary -p ${task.cpus} -e -n -v -z $gff
    """
}

Thanks in advance. Also, apologies as I can see why this might not strictly be bioinformatics, but it's where Nextflow has the most traction.

Workflow management has a firm foothold in bioinformatics and is in my opinion directly on topic on this site.

I have not read the full Nextflow manual in depth, but I think there is no Nextflow command to detect the maximum number of threads you can use.
However, you can try to automatize this task by using the nproc command.
nproc gives you the number of processing units available. An idea would be to include something such nproc - 1 as the max number of threads for Roary.

Thanks, as far I could see there are ways to dynamically allocate resources (https://www.nextflow.io/docs/latest/process.html#dynamic-computing-resources) but not detect threads, so this is a nice workaround.

Here's an example I wrote recently:
#!/usr/bin/env nextflow

params.requestedCPU = 80
params.requestedMem = 50

maxcpus = Runtime.runtime.availableProcessors()
requestedCpus = Math.round((params.requestedCPU * maxcpus) / 100)

maxmem = Runtime.runtime.totalMemory() / 10241024
requestedMem = (Math.round((params.requestedMem * maxmem) / 100))

process foo {
    echo true
    cpus = requestedCpus
    memory = requestedMem.GB

    script:
        """
        echo $requestedCpus
 echo $requestedMem
        """
}

You can pass the number of CPUs or Memory directly on the command line in percentage, e.g. nextflow run main.nf --requestedCPU 100 to have all CPU usages. Hope this helps.


text:
How to compute LD for pairs of variants with Plink
I have a set of SNPs (single nucleotide polymorphisms) { S1, S2, ..., SN } from approximately 200 humans. I wish to determine the linkage disequilibrium (LD) for a defined subset of these SNPs. By this I mean e.g. I would like to compute e.g. LD between S1 and S2, S3, S4, but not the LD between S2 and S3.
Could this calculation be performed in PLINK?

Hi Learner, thanks for you post and hope you can hang around here.The obvious question is why do you wish to sub-sample the genome, i.e. biological rationale?

In this case, I have a priori generated SNP pairs as the output of an eQTL analysis. Now I just want to know about the LD between pairs of SNPs, which I could compute manually, but I have thousands of pairs, and so it is ideal if there is efficient software like PLINK that I could use...

I've added a comment below

It could be that your question is not clear, but reads that you are wanting some SNPs subject to LD, but others excluded. For an LD to work you'll need a large number of SNPs which would be input simultaneously. You can subsample a given (human) population, your data set might therefore be further broken down from 200 humans.

Yes, you can do this with plink.
You can use the flag --ld-snp <variant_id> to calculate LD between your SNP of interest and all other SNPs.
If you have multiple SNPs of interest, and want each of those compared with every SNP in your dataset, use --ld-snps <snp1-snp100,snp102,snp104>
A final option to calculate LD between a specific pair of SNPs is --ld <variant_id1> <variant_id2>
This the relevant plink manual page

Anyway, it would be fairly easy via the filtering options here
The --exclude-snp  is what you are wanting.

I have worked in LD and subsampling a genome is open to forcing a result. The alternative approach is the subsample is the population, or examine disequilibrium with respect to the physical chromosome. A local shift in disequilibrium  or heterozygosity  (sliding window along a chromosome) is an indication of selective pressure. This approach is frequently used in eukaryotic parasites to identify drug resistance. I don't have experience of human genetics but there are certainly enough here who do. However, I think you need to define the biological rationale for the subsample. Just having loads of data is easily resolved because any algorithm can be written to handle automation.

The issue is that LD usually involves a 5% cut-off, even using a maximum likelihood approach I'd be cautious (I don't which PLINK uses), and is therefore sensitive variations in the data set


text:
How GFF3 attributes (9th column) varies from one gene prediction algorithm to another
GFF3 files are in tabular format with 9 fields per line, separated by tabs. The first 8 fields share almost same data structure, but the 9th field varies a lot depending on feature type and gene prediction algorithm.
Presently I am trying to build my own GFF module that will parse a GFF file and return one structured data object (a python dictionary structure). For that, I need a good understanding of the 9th field of the GFF format.
So can any one please help me by providing various GFF files?

You might want to look into this python GFF3 parser, which follows biopython's style: https://github.com/chapmanb/bcbb/tree/master/gff

Have you looked into gffutils? It will probably do most of the heavy lifting for you, i.e. parsing.

The first place to start is the GFF3 specification. This is the official word on what is and is not allowed in a GFF3 file. For example, users can define arbitrary attribute keys, so long as they do not begin with an uppercase letter (these are reserved for "official" use).
But your question doesn't seem to be about what is allowed, but what is commonly used. I have a question for you: which gene predictors are YOU using? Or are you using gene annotations produced by others?
Here are a few examples.
NCBI RefSeq
If you're using GFF3 from NCBI, chances are it looks something like this.
NC_007070.3     Gnomon  gene    329235  331223  .       +       .       ID=gene14;Dbxref=BEEBASE:GB42168,GeneID:551678;Name=LOC551678;gbkey=Gene;gene=LOC551678;gene_biotype=protein_coding
NC_007070.3     Gnomon  mRNA    329235  331223  .       +       .       ID=rna27;Parent=gene14;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;Name=XM_624067.4;gbkey=mRNA;gene=LOC551678;model_evidence=Supporting evidence includes similarity to: 66 ESTs%2C 24 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 113 samples with support for all annotated introns;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329235  329459  .       +       .       ID=id117;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329850  330082  .       +       .       ID=id118;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330166  330301  .       +       .       ID=id119;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330376  331223  .       +       .       ID=id120;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  CDS     329333  329459  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     329850  330082  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330166  330301  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330376  330416  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1

MAKER
The MAKER annotation workflow (paper, software) is a pretty commonly used gene annotation tool, and produces GFF3 output like this.
scaffold_12     maker   gene    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959
scaffold_12     maker   mRNA    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1;Parent=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959-mRNA-1;_AED=0.24;_eAED=0.18;_QI=0|0|0|0.66|0.5|0.33|3|0|218
scaffold_12     maker   exon    652527  652817  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1203;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    654877  655170  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1204;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    655272  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1205;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     652527  652817  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     654877  655170  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     655272  655343  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   gene    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851
scaffold_12     maker   mRNA    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1;Parent=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851-mRNA-1;_AED=0.22;_eAED=0.22;_QI=0|0|0.25|0.25|1|1|4|661|95
scaffold_12     maker   exon    941547  941631  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1206;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942343  942367  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1207;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942780  942920  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1208;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    943200  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1209;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     941547  941631  .       +       0       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942343  942367  .       +       2       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942780  942920  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     943200  943236  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   three_prime_UTR 943237  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:three_prime_utr;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1

GeneMark, SNAP, Augustus
If we look at the output of several common ab initio gene prediction tools, none of them actually produces GFF3 by default. Here is some output from GeneMark (paper, software).
543   1   +  Initial         799051     799097         47          1 2 - -
543   2   +  Terminal        799236     799266         31          3 3 - -

544   3   -  Terminal        802357     802514        158          3 2 - -
544   2   -  Internal        802607     802685         79          1 1 - -
544   1   -  Initial         802829     802843         15          3 1 - -

Here is some output from SNAP (paper, software).
Einit   8230    8239    +       9.329   0       1       0       scaffold_12-snap.4
Exon    8848    8869    +       2.484   2       2       2       scaffold_12-snap.4
Exon    10121   10208   +       15.302  1       0       2       scaffold_12-snap.4
Exon    11361   11420   +       5.969   0       0       2       scaffold_12-snap.4
Exon    11471   11535   +       0.921   0       2       1       scaffold_12-snap.4
Eterm   12169   12187   +       18.163  1       0       1       scaffold_12-snap.4
Einit   14569   14668   +       -6.918  0       1       0       scaffold_12-snap.5
Exon    15029   15203   +       -1.023  2       2       0       scaffold_12-snap.5
Exon    16171   16307   +       -9.230  1       1       1       scaffold_12-snap.5
Eterm   16667   16698   +       9.829   2       0       0       scaffold_12-snap.5
Einit   17809   17898   -       -5.390  0       0       0       scaffold_12-snap.6
Exon    17299   17350   -       6.978   0       1       2       scaffold_12-snap.6
Eterm   17029   17129   -       -2.857  2       0       0       scaffold_12-snap.6

Both of these formats appear to be tab-delimited, but that's where the similarity with GFF3 ends.
The output of Augustus (paper, software) looks like GTF (a variant of the GFF format) once you lose all the lines beginning with a # symbol, but Augustus also has a GFF3 output mode.
# Predicted genes for sequence number 1 on both strands
# start gene g1
scaffold_12     AUGUSTUS        gene    2841    3400    1       +       .       g1
scaffold_12     AUGUSTUS        transcript      2841    3400    1       +       .       g1.t1
scaffold_12     AUGUSTUS        start_codon     2841    2843    .       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        intron  3027    3097    1       +       .       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     2841    3026    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     3098    3400    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        stop_codon      3398    3400    .       +       0       transcript_id "g1.t1"; gene_id "g1";
# protein sequence = [MAIKNAEHDLRVIVDAIEGLGLKVAPHKTEAMAFPASALCGRRGAAPPKIRLGGSSILVGSRSRWYISGHSENSSKSP
# RTEGKETTPLQQRDPLDAPLWVSGVVAHCCGGPEGQEGCPGLAAQGSDQGVLRIRDGLLCGYDGCGDHRPRPSDSSAGGGLCRP]
# end gene g1
###
# start gene g2
scaffold_12     AUGUSTUS        gene    4712    15229   0.21    +       .       g2
scaffold_12     AUGUSTUS        transcript      4712    15229   0.21    +       .       g2.t1
scaffold_12     AUGUSTUS        start_codon     4712    4714    .       +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  4858    5591    0.74    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  5686    15028   0.23    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     4712    4857    0.72    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     5592    5685    0.42    +       1       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     15029   15229   0.43    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        stop_codon      15227   15229   .       +       0       transcript_id "g2.t1"; gene_id "g2";
# protein sequence = [MGRNSHRSCCVVNCKITSAKSDCKFYKFPTAKWKINQRKMWVAAVKRQKYIKDEISHAETQTEITEVTGATKVNYANK
# KYICLLFVRTYVLRMLVDVALSNLRFSLFGIRKSLEIFGQSEKADQTRWRLPSCEMEWIESRKGKMRE]
# end gene g2
###

Parsing attributes
So...how should YOU handle attributes in GFF3's 9th column? That depends a lot on what you want to do with the data.
The most important attributes are ID and Parent, which are used to define relationships between features and subfeatures. (These relationships implicitly define a directed acyclic graph of features, although most GFF3 parsers don't directly support traversal of this graph.)
But not only do you have to handle the pre-defined attributes discussed in the GFF3 specification, you also have to be able to handle any number of arbitrary attributes whose keys you may not know beforehand.
The simplest way to handle this would be to parse the attribute column into a dictionary of key/value pairs. Once it's in this form, it's trivial to see what attributes are there and how to access them.
For example, if we have a feature with the following attributes...
ID=mRNA42;Parent=gene19;integrity=0.95;foo=bar

...we would want it in a dictionary like so.
attributes = {
    'ID': 'mrna42',
    'Parent': 'gene19',
    'integrity': '0.95',
    'foo': 'bar',
}

Python code to parse that might look something like this.
attributes = dict()
for keyvaluepair in attributestring.split(';'):
    for key, value in keyvaluepair.split('='):
        attributes[key] = value

One thing to consider is that it's possible for an attribute to have multiple values (separated by commas). It's not commonly used, but it is valid and something that should be handled or at least checked for.

Thank's @Daniel for your excellent answer.My ultimate goal is to convert any gff3 format to a python dictionary structure. can you please help me convert one gene record(One gene record may have mRNA, exon, CDS etc) into a dictionary structure?

@Arijit Parsing a single line into a dictionary is trivial, but resolving parent/child relationships involves ID cross-referencing and is more complicated. Most GFF3 parsers leave this as an exercise to the coder. I'm working on a Python library called tag right now. Feel free to check it out and see if it meets your needs. http://tag.readthedocs.io/


text:
Removing PCR duplicates in RNA-seq Analysis
After reading some of the forum posts in Biostar and SeqAnswers I find it very confusing whether to filter out the duplicate reads from aligned files or not. As far I understand it's very difficult to distinguish between highly expressed genes and duplicate reads and we may lose important information during the filtration process.
So, is it really necessary to remove the duplicates in differential expression analysis using RNA-seq data?

It all depends on how severe the problem is. In general you should only remove them if it is a very high percentage.

Any particular threshold like more than 35%?

Unfortunatly not. The problem is by removing them you might remove real duplicated data thereby under powering your analysis - but by not removing them you might count PCR artefacts thereby increasing false positive. So either option have a potential bias.
I would probably try doing your analysis both with and without removal and then manually inspect some of the results only found by one approach to se which you find more trustworthy. What you should pay special attention to are whether the coverage is more uniform or seem to be affected by loads of reads from single positions.

For normal RNA-seq PCR duplicates are normally kept in, but the duplication rate can be used as a quality control: The higher the duplication rate, the lower the quality. For expression analysis, it is probably best to discard high duplication rate samples, rather than deduplicate them.
In general, the smaller the amount of RNA input into the library prepartion the worst the duplication. Many protocols for very low input quantities (such as single cell) include random barcodes called UMIs (Unique Molecular Identifiers). These allow PCR duplicates to be distinguished from genuinely independent molecules that just happen to have the sample mapping position.

Since the OP does not mention which RNA classes the analysis, let me add that UMIs are also extremely important for sRNA, in particular piRNAs, which quite often map to the exact same locations in the genome and have the same sequence. As afar I know current methods to mark duplicates in mapped reads would result in extreme under-reporting of these RNAs.

Generally you should just leave them as is. One does remove/mark duplicates in DNA seq.
For further read check this Nature paper

Your answer could be improved by hinting on why, according to the paper you cite, removing "PCR duplicates" shouldn't be done in RNASeq. E.g., by citing the one crucial sentence of the abstract: "We find that a large fraction of computationally identified read duplicates are not PCR duplicates and can be explained by sampling and fragmentation bias."

Sure..I could have excuse me for short answer.

Some of the researchers in my institute have researched on this and use UMIs (Unique Molecular Identifiers) to eliminate PCR duplicates. They have been able to increase the reproducibility of the rna seq analysis they had been working on. The UMIs capture some lengths of diverse nucleotides from the rna-sequence. I cannot validate how much UMIs can help in deduplication, as even Nature journal says that UMIs specifically work on the 3’ end. But it has been observed that eliminating duplicates doesn’t improve the accuracy of quantification and it could introduce a bias in the data. So do not eliminate the duplicates, it's just a waste of time.

Hello Shirin, I think it's important to provide a reference for the following statement in your answer: "But it has been observed that eliminating duplicates doesn’t improve the accuracy of quantification and it could introduce a bias in the data."


text:
Is there public RESTful api for Gnomad?
I currently find Harvard's RESTful API for ExAC extremely useful and I was hoping that a similar resource is available for Gnomad?
Does anyone know of a public access API for Gnomad or possibly any plans to integrate Gnomad into the Harvard API?

Just to add a comment that Cellbase annotation now includes Gnomad data (Exomes and Genomes) and can be accessed via RESTful API -
http://bioinfo.hpc.cam.ac.uk/cellbase/webservices

As far as I know, no but the vcf.gz files are behind a http server that supports Byte-Range, so you can use tabix or any related API:
$ tabix "https://storage.googleapis.com/gnomad-public/release-170228/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz" "22:17265182-17265182"
22  17265182    .   A   T   762.04  PASS    AC=1;AF=4.78057e-06;AN=209180;BaseQRankSum=-4.59400e+00;ClippingRankSum=2.18000e+00;DP=4906893;FS=1.00270e+01;InbreedingCoeff=4.40000e-03;MQ=3.15200e+01;MQRankSum=1.40000e+00;QD=1.31400e+01;ReadPosRankSum=2.23000e-01;SOR=9.90000e-02;VQSLOD=-5.12800e+00;VQSR_culprit=MQ;GQ_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1;DP_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0|0;AB_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0;GQ_HIST_ALL=1591|589|120|301|650|589|1854|2745|1815|4297|5061|2921|10164|1008|6489|1560|7017|457|6143|52950;DP_HIST_ALL=2249|1418|6081|11707|16538|9514|28624|23829|7391|853|95|19|1|0|0|1|0|1|0|0;AB_HIST_ALL=0|0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0;AC_AFR=0;AC_AMR=0;AC_ASJ=0;AC_EAS=0;AC_FIN=1;AC_NFE=0;AC_OTH=0;AC_SAS=0;AC_Male=1;AC_Female=0;AN_AFR=11994;AN_AMR=31324;AN_ASJ=7806;AN_EAS=13112;AN_FIN=20076;AN_NFE=94516;AN_OTH=4656;AN_SAS=25696;AN_Male=114366;AN_Female=94814;AF_AFR=0.00000e+00;AF_AMR=0.00000e+00;AF_ASJ=0.00000e+00;AF_EAS=0.00000e+00;AF_FIN=4.98107e-05;AF_NFE=0.00000e+00;AF_OTH=0.00000e+00;AF_SAS=0.00000e+00;AF_Male=8.74386e-06;AF_Female=0.00000e+00;GC_AFR=5997,0,0;GC_AMR=15662,0,0;GC_ASJ=3903,0,0;GC_EAS=6556,0,0;GC_FIN=10037,1,0;GC_NFE=47258,0,0;GC_OTH=2328,0,0;GC_SAS=12848,0,0;GC_Male=57182,1,0;GC_Female=47407,0,0;AC_raw=1;AN_raw=216642;AF_raw=4.61591e-06;GC_raw=108320,1,0;GC=104589,1,0;Hom_AFR=0;Hom_AMR=0;Hom_ASJ=0;Hom_EAS=0;Hom_FIN=0;Hom_NFE=0;Hom_OTH=0;Hom_SAS=0;Hom_Male=0;Hom_Female=0;Hom_raw=0;Hom=0;POPMAX=FIN;AC_POPMAX=1;AN_POPMAX=20076;AF_POPMAX=4.98107e-05;DP_MEDIAN=58;DREF_MEDIAN=5.01187e-84;GQ_MEDIAN=99;AB_MEDIAN=6.03448e-01;AS_RF=9.18451e-01;AS_FilterStatus=PASS;CSQ=T|missense_variant|MODERATE|XKR3|ENSG00000172967|Transcript|ENST00000331428|protein_coding|4/4||ENST00000331428.5:c.707T>A|ENSP00000331704.5:p.Phe236Tyr|810|707|236|F/Y|tTc/tAc||1||-1||SNV|1|HGNC|28778|YES|||CCDS42975.1|ENSP00000331704|Q5GH77||UPI000013EFAE||deleterious(0)|benign(0.055)|hmmpanther:PTHR14297&hmmpanther:PTHR14297:SF7&Pfam_domain:PF09815||||||||||||||||||||||||||||||,T|regulatory_region_variant|MODIFIER|||RegulatoryFeature|ENSR00000672806|TF_binding_site|||||||||||1||||SNV|1||||||||||||||||||||||||||||||||||||||||||||,T|regulatory_region_variant|MODIFIER|||RegulatoryFeature|ENSR00001729562|CTCF_binding_site|||||||||||1||||SNV|1||||||||||||||||||||||||||||||||||||||||||||

UPDATE: 2019: the current server for gnomad doesn't support Byte-Range requests.

As of 15 Jan 2020, cellbase (bioinfo.hpc.cam.ac.uk/cellbase/webservices) does. See Pasted's comment on the question. Byte-range requests are supported, and results come from GnomAD, ExAC, and other data sets. A pretty awesome resource actually, even if byte-range requests were supported at GnomAD or UCSC.

The new gnomAD site (as of August 2019) says no, no API yet:

How do I query a batch of variants? Do you have an API?
We currently do not have a way to submit batch queries on the browser, but we are actively working on developing an API for ExAC/gnomAD. If you would like to learn about GraphQL, which we will use to work with the API, an overview can be found at https://graphql.org.  You can also obtain information on all variants from the VCFs and Hail Tables available on our downloads page.

But, the web interface itself already makes POST requests to https://gnomad.broadinstitute.org/api to send and receive JSON/GraphQL.  So, you can make those same queries programmatically right now, even if it's not officially a public API.
Here's an example in Python to get some basic info on variants for a particular gene.  This way you get simple nested Python objects to work with:
  { 'consequence': 'intron_variant',
    'pos': 6928442,
    'rsid': 'rs782435448',
    'variant_id': '12-6928442-C-A'},
  { 'consequence': 'splice_region_variant',
    'pos': 6928462,
    'rsid': None,
    'variant_id': '12-6928462-C-A'},
  { 'consequence': 'splice_acceptor_variant',
    'pos': 6928464,
    'rsid': 'rs782577109',
    'variant_id': '12-6928464-G-A'},
  { 'consequence': 'missense_variant',
    'pos': 6928466,
    'rsid': 'rs782208003',
    'variant_id': '12-6928466-C-T'},

(I found it useful to go this route because then the full metadata visible in the gnomAD web interface is then available, including the per-variant details like allele counts by population.  I couldn't find this information in the other APIs described here.)

can you expand on your example a bit. the graphAPI is unfamiliar to me. How would you turn the json above into an actual query?

That JSON is the output for the GraphQL used in the linked example code.  (Things may have changed since then, too, since that was back in 2019 and I haven't dealt with gnomAD since.)

You can browse gnomAD variants with ClinGen Allele Registry (there is API spec available).

I faced same issue recently, I found those link and python script:
gnomAD GraphQL api
https://gnomad.broadinstitute.org/api
It works great but it is a kind of different query language.
Please check here for the docs: https://graphql.org/learn/queries/
gnomAD Python Api
https://github.com/furkanmtorun/gnomad_python_api

I found Jesse's code quite usefull ! For those who try to reproduce it, you should now add the reference genome ID, such as :
#!/usr/bin/env python

import requests
# import pprint
# prettyprint = pprint.PrettyPrinter(indent=2).pprint

def fetch(jsondata, url="https://gnomad.broadinstitute.org/api"):
    # The server gives a generic error message if the content type isn't
    # explicitly set
    headers = {"Content-Type": "application/json"}
    response = requests.post(url, json=jsondata, headers=headers)
    json = response.json()
    if "errors" in json:
        raise Exception(str(json["errors"]))
    return json

def get_variant_list(gene_id, dataset="gnomad_r2_1"):
    # Note that this is GraphQL, not JSON.
    fmt_graphql = """
    {
        gene(gene_id: "%s", reference_genome: GRCh38) {
          variants(dataset: %s) {
            consequence
            pos
            rsid
            variant_id: variantId
          }
        }
      }
    """
    # This part will be JSON encoded, but with the GraphQL part left as a
    # glob of text.
    req_variantlist = {
        "query": fmt_graphql % (gene_id, dataset),
        "variables": {"withFriends": False}
        }
    response = fetch(req_variantlist)
    return response["data"]["gene"]["variants"]

prettyprint(get_variant_list("ENSG00000010610"))

Thanks, I updated my code to include the reference genome as well (you're right, it no longer works without that).

I created a python package based on SQLite databases, where you can easily query all gnomAD variants for GRCh37/38. https://github.com/KalinNonchev/gnomAD_DB I have precomputed SQLite databases for gnomAD WGS for GRCh37/38 in the description of the package. Please take a look there."

https://github.com/KalinNonchev/gnomAD_MAF would solve your problem I think.
It is in python and you can annotate fast variants with their allele frequencies.


text:
The state, limitations and comparisons of large variant stores
Background: We're increasingly needing some way of storing lots of variant data associated with lots of subjects: think clinical trials and hospital patients, looking for disease-causing or relevant genes. A thousand subjects is where we'd start, there's talk of millions on the horizon. With various genomic medicine initiatives, this is likely a wider need.
The problem: While there's plenty of platforms out there, it's a rapidly evolving field. It's difficult to get a feel for how (and if) they perform and how they line up against each other:

What's scalable and can handle a lot of data? What sort of limits?
What's robust and not a teetering pile of hacked-together components?
What has a large community behind it and is actually used widely?
What makes for easy access and search from another service? (Commandline, REST or software APIs)
What sort of variants they handle?
What sort of parameters can be used in searching?

Solutions I've seen so far:

BigQ: used with i2b2, but its wider use is unclear
OpenCGA: looks the most developed, but I've heard complaints about the size of data it spits out
Using BigQuery over a Google Genomics db: doesn't seem to be a general solution
Gemini: recommended but is it really scalable and accessible from other services?
SciDb: a commercial general db
Quince
LOVD
Adam
Whatever platform DIVAS & RVD run on: which may not be freely available
Several graphical / graph genome solutions: We (and most other people) are probably not dealing with graph genome data at the moment, but is this a possible solution?
Roll your own: Frequently recommended but I'm sceptical this is a plausible solution for a large dataset.

Anyone with experience give a review or high-level guide to this platform space?

My two cents: use MongoDB wrapped in a simple REST framework.  Allows flexible model and queries and should scale into the billions of records on a single node.  Working on a FLOSS project for this at the moment, but it is not production ready yet.

@woemler How is it compared to other approaches? Someone I know tried MongoDB ~5 years ago on 1000g genotypes. He said MongoDB was over 10x slower than bcf2 on parallel queries while having a much larger disk/memory footprint. That said, he was new to MongoDB back then and might not be doing it the optimal way.

@user172818: The newer versions of MongoDB (3.2+) are significantly faster than the versions from several years ago.  I have benchmarked it against other free RDBMS's, and it typically performs as-well-as or better, especially for complex data representations, like variant calls.

Is storing the data more important here, or is processing statistics (using Python, R, etc..) about the data more important?

@macgyver : good observation. The data - it's supposedly people will want to be mining and querying the data, rather than looking at summary stats and analyses.

An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.
For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in addition to genotypes. Hail is at least something heavily used in practice, although mostly by a few groups so far.
A simpler problem is to store genotypes only. This is sufficient to the majority of end users. There are better approaches to store and query genotypes. GQT, developed by the Gemini team, enables fast query of samples. It allows you to quickly pull samples under certain genotype configurations. As I remember, GQT is orders of magnitude faster than google genomics API to do PCA. Another tool is BGT. It produces a much smaller file and provides fast and convenient queries over sites. Its paper talks about ~32k whole-genome samples. I am in the camp who believe specialized binary formats like GQT and BGT are faster than solutions built on top of generic databases. I would encourage you to have a look if you only want to query genotypes.
Intel's GenomicDB approaches the problem in a different angle. It does not actually keep a "squared" multi-sample VCF internally. It instead keeps per-sample genotypes/annotations and generates merged VCF on the fly (this is my understanding, which could be wrong). I don't have first-hand experience with GenomicDB, but I think something in this line should be the ultimate solution in the era of 1M samples. I know GATK4 is using it at some step.
As to others in your list, Gemini might not scale that well, I guess. It is partly the reason why they work on GQT. Last time I checked, BigQuery did not query individual genotypes. It only queries over site statistics. Google genomics APIs access individual genotypes, but I doubt it can be performant. Adam is worth trying. I have not tried, though.

+1 for hail, clearly The Right Answer at this point

You can query individual genotypes using BigQuery.  The biggest challenge at this point is having to write your own queries to do analysis.


text:
BAM to BigWig without intermediary BedGraph
I have a pipeline for generating a BigWig file from a BAM file:
BAM -> BedGraph -> BigWig

Which uses bedtools genomecov for the BAM -> BedGraph part and bedGraphToBigWig for the BedGraph -> BigWig part.
The use of bedGraphToBigWig to create the BigWig file requires a BedGraph file to reside on disk in uncompressed form as it performs seeks. This is problematic for large genomes and variable coverage BAM files when there are more step changes/lines in the BedGraph file. My BedGraph files are in the order of 50 Gbytes in size and all that IO for 10-20 BAM files seems unnecessary.
Are there any tools capable generate BigWig without having to use an uncompressed BedGraph file on disk? I'd like this conversion to hapen as quickly as possible.
I have tried the following tools, but they still create/use a BedGraph intermediary file:

deepTools

Some Benchmarks, Ignoring IO
Here are some timings I get for creating a BigWig file from a BAM file using 3 different pipelines. All files reside on a tmpfs i.e. in memory.
BEDTools and Kent Utils
This is the approach taken by most.
time $(bedtools genomecov -bg -ibam test.bam -split -scale 1.0 > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    1m20.015s
user    0m56.608s
sys     0m27.271s

SAMtools and Kent Utils
Replacing bedtools genomecov with samtools depth and a custom awk script (depth2bedgraph.awk) to output bedgraph format has a significant performance improvement:
time $(samtools depth -Q 1 --reference test.fasta test.bam \
  | mawk -f depth2bedgraph.awk \
  > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    0m28.765s
user    0m44.999s
sys     0m1.166s

Although it is has less features, we used mawk here as it's faster than gawk (we don't need those extra features here).
Parallelising with xargs
If you want to have a BigWig file per chromosome, you can easily parallelise this across chromosome/reference sequences. We can use xargs to run 5 parallel BAM->BedGraph->BigWig pipelines, each using the tmpfs mounted /dev/shm for the intermediary BedGraph files.
cut -f1 test.fasta.chrom.sizes \
  | xargs -I{} -P 5  bash -c 'mkdir /dev/shm/${1} \
  && samtools depth -Q 1 --reference test.fasta -r "${1}" test.bam \
  | mawk -f scripts/depth2bedgraph.awk \
  > "/dev/shm/${1}/test.bam.bedgraph" \
  && mkdir "./${1}" \
  && bedGraphToBigWig \
    "/dev/shm/${1}/test.bam.bedgraph" \
    test.fasta.chrom.sizes \
    "./${1}/test.bam.bw" \
  && rm "/dev/shm/${1}/test.bam.bedgraph"' -- {}

deepTools
Let's see how deepTools performs.
time bamCoverage --numberOfProcessors max \
  --minMappingQuality 1 \
  --bam test.bam --binSize 1 --skipNonCoveredRegions \
  --outFileName deeptools.bw

real    0m40.077s
user    3m56.032s
sys     0m9.276s

It's been on my to do list to remove the "make a single bedGraph file" intermediate step from deepTools. Given how that's all implemented there's no good way to avoid at least some intermediate file creation (without blowing up the memory at least).

I needed to generate the intermediate bedgraph files when I was creating BigWig files for GBrowse in 2011; it's disappointing to see that it hasn't changed much in the many years since then. The file seeking is a killer -- it excludes most of the shell tricks that can be used to work around programs that don't explicitly support the common UNIX pipeline (e.g. <(), FIFO buffers, /dev/stdin).

@DevonRyan when you say blow up memory, will that be equivilent to the size of the intermediary bedgraph file? I have a 720Gbyte fat memory compute node...so for me, speed is a bigger issue than memory. Perhaps an option to state the max memory usage, at which point you start going to disk?

@NathanS.Watson-Haigh Yeah, it'd be equivalent more or less to the file size. deepTools itself was never engineered to make this easier to implement and the generally crappy nature of multithreading in python doesn't help that :(

@DevonRyan how about using the system's temp dir as the default location for the intermediary bedgraph file and expose this as a command line argument. That way I could use my raided SSD's for more convenient and faster IO. I was doing this with my bedtools and bedGraphToBigWig pipeline.

The system temp dir (or shared memory) is already what's used for all of the tmp files. Honestly, the current implementation is a hold-over from before I wrote pyBigWig (or even joined the group) and deepTools was using bedGraphToBigwig in a system call. It's mostly a matter of me finding time to change how bamCoverage does this final step.

@DevonRyan is there any reason a named pipe as suggested by my answer wouldn't work for this? I have never used these tools so I can't be sure.

@terdon: Sorry, I was referring specifically to deepTools, which was mentioned in the post. One can certainly do what you've described, it'll just be slower than what deeptools can do (but better with IO, so maybe it's all a wash at that point!).

@terdon I don't see your named pipe based answer

@bli no I deleted it since it won't work here because the question requires the ability to seek the file. FIFOs won't work with seeking AFAIK, so the answer wasn't useful. I hadn't realized the OP required a seekable file.

This can be done in R very easily from an indexed .bam file.
Given single-end file for sample1.
library(GenomicAlignments)
library(rtracklayer)

## read in BAM file (use readGAlignmentPairs for paired-end files)
gr <- readGAlignments('sample1.bam')

## convert to coverages
gr.cov <- coverage(gr)

## export as bigWig
export.bw(gr.cov,'sample1.bigwig')

Be aware that this method doesn't include normalization steps (such as normalizing to total coverage). Most of these additional steps can be added if necessary.

I had hoped for a command line tool but I'm not totally against going into R if needs be. I need to test more extensively with regards to speed so I need to have more of a play. If tests show that it's slow for the number and sizes of BAM files I have, or a better command line answer materialises, I will accept this answer.

Just for reference, there is a not very in depth, but first-hand reasoning (Jim Kent) given here about why bedGraphToBigWig does not support streaming http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin
Edit: as noted wigToBigWig allows bedgraph streaming input via stdin but takes (possibly much) more memory than bedGraphToBigWig itself does
cat file.bedgraph | wigToBigWig stdin chrom.sizes output.bigwig

Note: You can also use wiggletools bin 10 for example to bin rows and downsample

Ah! "However if you do have enough memory, you can still use wigToBigWig on bedGraphs, and it will fit in a pipe."

That's sort of a weird decision. The zoom levels are written after the raw entries, so one can write all of the raw interval entries to disk and then do another pass over them when making the zoom levels. Then you don't need to store everything in memory. That's how libBigWig works, in fact.

The QuasR package in R has a useful function called qExportWig.
To use it you first need to create a sample sheet and a qProject object - this can be done starting with bam files or fastq files (in which case it also  performs alignments). You can then create wig or bigWig files without writing a bedGraph file to disk.
An example would look something like this:
library(QuasR)

# Create qProject object using qAlign
# Sample file should contain tab-separated file names (BAM or fastq) and sample names
myproj <- qAlign("samples.txt", genome=mygenome, paired="no")

# Export all samples in project as bigWig, scaling to mean number of alignments by default
qExportWig(myproj, binsize=50, createBigWig=T)

It has a lot of flexible options, including scaling or log-transforming counts, shifting reads and changing the bin size.

It would be helpful to provide some example code to demonstrate how one can go from BAM to BigWig. Like the answer provided by @story


text:
Transcript Coordinate Ranges to Genomic Coordinates
I have 2 GFF3 files:

Features using transcript IDs as the landmarks. i.e. "CDS" feature types using coordinates from transcript space.
Features using chromosome IDs as the landmarks. i.e. "exon" feature types using coordinates from chromosome space.

I would like to transform the coordinates of the features in file 1 into the coordinate space of file 2. i.e. transcript-based coordinates into genomic-based coodinates.
Here is an example for file1:
cat transcript_orfs.gff3
##gff-version 3
##sequence-region Tx.1 1 4000
Tx.1    ORF_finder  gene    1   4000    .   +   .   ID=1
Tx.1    ORF_finder  mRNA    1   4000    .   +   .   ID=2;Parent=1
Tx.1    ORF_finder  exon    1501    2500    .   +   0   ID=3;Parent=2
Tx.1    ORF_finder  CDS 1501    2500    .   +   0   ID=4;Parent=2

gt sketch -addintrons transcript_orfs.png transcript_orfs.gff3

Here is an example for file1:
cat genome.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 genome_assembly gene    1001    6000    .   +   .   ID=Tx
chr3A_part1 genome_assembly mRNA    1001    6000    .   +   .   ID=Tx.1;Parent=Tx
chr3A_part1 genome_assembly exon    1001    3000    .   +   .   ID=Tx.1.exon1;Parent=Tx.1
chr3A_part1 genome_assembly five_prime_UTR  1001    2000    .   +   .   ID=Tx.1.utr5;Parent=Tx.1
chr3A_part1 genome_assembly CDS 2001    3000    .   +   0   ID=Tx.1.cds1;Parent=Tx.1
chr3A_part1 genome_assembly exon    4001    6000    .   +   .   ID=Tx.1.exon2;Parent=Tx.1
chr3A_part1 genome_assembly CDS 4001    5000    .   +   2   ID=Tx.1.cds2;Parent=Tx.1
chr3A_part1 genome_assembly three_prime_UTR 5001    6000    .   +   .   ID=Tx.1.utr3;Parent=Tx.1

gt sketch -addintrons genome.png genome.gff3

I want to convert the coordinates of the features from file1 into chromosome genomic coordinates. I expect to get something like the following output:
cat output.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 ORF_finder  gene    1001    6000    .   +   .   ID=1
chr3A_part1 ORF_finder  mRNA    1001    6000    .   +   .   ID=2;Parent=1
chr3A_part1 ORF_finder  exon    2501    3000    .   +   0   ID=3.1;Parent=2
chr3A_part1 ORF_finder  CDS 2501    3000    .   +   0   ID=4.1;Parent=2
chr3A_part1 ORF_finder  exon    4001    4500    .   +   0   ID=3.2;Parent=2
chr3A_part1 ORF_finder  CDS 4001    4500    .   +   1   ID=4.2;Parent=2

gt sketch -addintrons output.png output.gff3

I have looked into using mapFromTranscripts() from Bioconductor's GenomicRanges library but I have made little progress trying to decipher the manual.

This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances.
R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're going. :-)
#!/usr/bin/env python3

from argparse import ArgumentParser, FileType
from collections import defaultdict
from re import search, sub

cli = ArgumentParser()
cli.add_argument('txspace', type=FileType('r'))
cli.add_argument('chrspace', type=FileType('r'))
args = cli.parse_args()

exons = defaultdict(list)
for line in args.txspace:
    if line.startswith('#') or line.strip() == '':
        continue
    values = line.split('\t')
    featuretype = values[2]
    if featuretype not in ('CDS', 'exon'):
        continue
    txid = values[0]
    exons[txid].append(values)

for line in args.chrspace:
    values = line.split('\t')
    if len(values) != 9:
        print(line, end='')
        continue
    featuretype = values[2]
    if featuretype not in ('gene', 'mRNA'):
        continue
    print(line, end='')
    if featuretype == 'mRNA':
        txid = search('ID=([^;\n]+)', values[8]).group(1)
        start = int(values[3])
        txexons = exons[txid]
        for exon in txexons:
            exon[0] = values[0]
            exon[3] = str(int(exon[3]) + start)
            exon[4] = str(int(exon[4]) + start)
            exon[8] = sub(r'Parent=[^;\n]+', 'Parent=' + txid, exon[8])
            print('\t'.join(exon), end='')

The linked gist contains a class I wrote for doing this on GTF files:
https://gist.github.com/IanSudery/d8349c22823a475ceb489c3e8aeb448e
It uses the GTF class from cgat, which can be found here:
https://github.com/cgat-developers/cgat-apps
You would use it to do this task like so:
from cgat import gtf
from cgat import iotools
import TranscriptCoordsInterconverter

trans_coord_file = gtf.iterator(iotools.open_file(args[1])
transcripts = dict()

for transcript in gtf.transcript_iterator(trans_coords_file):
    transcript_id = transcript[0].transcript_id
    transcripts[transcript_id] = transcript

trans_coord_file.close()

genome_coord_file = gtf.iterator(iotools.open_file(args[2]))
for transcript in gtf.transcript_iterator(genome_coord_file):

    transcript_id = transcript[0].transcript_id
    contig = transcript[0].contig
    strand = transcript[0].strand

    if transcript_id not in transcripts:
        continue

    converter = TranscriptCoordsInterconverter(transcript)
    for exon in transcripts[transcript_id]:
        genome_intervals = converter.transcript_interval2genome_intervals((exon.start,exon.end))
        for interval in genome_intervals:
            new_entry = GTF.Entry().fromGTF(exon)
            new_entry.contig = contig
            new_entry.strand = strand
            new_entry.start = interval[0]
            new_entry.end = interval[1]
            print (str(new_entry)+"\n")