Question

makeTxDbFromGFF drops genes which have multiple chromosome locations. (with iGenome GTF)

1

Entering edit mode

Marlin ▴ 20

@marlin-11371

Last seen 6.3 years ago

For a GTF file downloaded from iGenome， say, hg19.gtf, a significant proportion of genes will be dropped without warning when use makeTxDbFromGFF to import it.

As I examined those genes in the source GTF file, I believe the genes with multiple chromosome locations (e.g. chr1 and chr2) will be dropped when using makeTxDbFromGFF.

For example, if you grep a gene name called "XGY2" in the iGenome GTF file, you get:

chrX    unknown    exon    2670337    2670376    .    +    .    gene_id "XGY2"; gene_name "XGY2"; transcript_id "NR_003254"; tss_id "TSS24556";
chrX    unknown    exon    2688591    2688632    .    +    .    gene_id "XGY2"; gene_name "XGY2"; transcript_id "NR_003254"; tss_id "TSS24556";
chrX    unknown    exon    2692757    2693037    .    +    .    gene_id "XGY2"; gene_name "XGY2"; transcript_id "NR_003254"; tss_id "TSS24556";
chrY    unknown    exon    2620337    2620376    .    +    .    gene_id "XGY2"; gene_name "XGY2"; transcript_id "NR_003254_1"; tss_id "TSS3943";
chrY    unknown    exon    2638591    2638632    .    +    .    gene_id "XGY2"; gene_name "XGY2"; transcript_id "NR_003254_1"; tss_id "TSS3943";
chrY    unknown    exon    2642757    2643037    .    +    .    gene_id "XGY2"; gene_name "XGY2"; transcript_id "NR_003254_1"; tss_id "TSS3943";

However, this gene will not be imported into the TxDb using makeTxDbFromGFF, it is just silently dropped.

I don't know if anyone has faced the problem before. It could be the case that TxDb is designed using unique gene id, but I would consider it as a severe bug as people are losing their annotations without awareness.

Please help to reproduce the problem, thanks.

genomicfeatures gtf maketxdbfromgff iGenome • 2.4k views

ADD COMMENT • link updated 8.0 years ago by Martin Morgan 25k • written 8.0 years ago by Marlin ▴ 20

score 2 · Accepted Answer · 2016-11-13

Marlin provided a nice reproducible example on the Bioc-devel mailing list. The problem is not that genes on multiple chromosomes are not imported, but rather that genes() by default reports only genes that have all exons on the same strand. In the current GenomicFeatures (1.26.0) the help page ?genes shows

     genes(x, columns="gene_id", filter=NULL, vals=NULL,
              single.strand.genes.only=TRUE)
...
single.strand.genes.only: 'TRUE' or 'FALSE'. If 'TRUE' (the default),
          then genes that have exons located on both strands of the
          same chromosome or on two different chromosomes are dropped.
          In that case, the genes are returned in a GRanges object.
          Otherwise, all genes are returned in a GRangesList object
          with the columns specified thru the 'columns' argument set as
          _top level_ metadata columns. (Please keep in mind that the
          _top level_ metadata columns of a GRangesList object are not
          displayed by the 'show' method.)

I believe this is because the original contract of genes() was to return a GRanges() object, which is not consistent with exons on separate chromosomes. Other functions, e.g., exons(), did not face this constraint and report all exons.

A related recent change (in response to keepSeqlevels for GRangesList: drop elements with multiple seqnames ) adds pruning.mode to seqlevels<-, so that genes with exons on several seqlevels are not necessarily dropped; I believe this change is still being incorporated.