[sindrle]

Thank you, you have been very helpful.

Im currently struggling with how to implement featureCounts to DEXseq, its much easier to use HTSeq, since its described in the manual.

I really want to be able to use featureCounts, since with HTSeq I had to "hack" the python script to work with my annotation file, whereas I don't need that with featureCounts.

Also, I want to be able, in the end, to present expression values of transcripts as counts/library size/gene length, which also is easier with featureCounts, since it outputs gene length!

[yangliao]

The values in the "Gene Length" column are the length of each feature (on feature mode) or the total length of all the features belonging to each gene (on meta-feature mode).

It is straightforward on the feature mode: gene length = stop - start + 1. However, on the meta-feature mode, the total length is the uniquely covered length of the gene; namely, if two exons in this gene have an overlapping part, the overlapping part is measured only once in the Gene Length column.

For example, on the meta-feature mode, if a gene has three exons: [100, 200], [100, 300] and [200, 250], the Gene Length value for this gene is only 201 because the three exons overlap with each other, and they can be merged into one interval on the chromosome: [100, 300], which has 201 bases.

Ah! Now I get it, very nice.

One more question, where is the information about "gene length" stored?

[sindrle]

Ah! Now I get it, very nice.

One more question, where is the information about "gene length" stored?

[yangliao]

Hi Sindrle,

If you look at the annotation component in the returned R object, you will find another table containing the details of the annotations.

Each row in adiExons$counts corresponds to the same row in adiExons$annotation, so you can tell which count number is for which exon. The information in the annotation table (e.g., chromosome name, start, stop, gene name, strand symbol and exon length) can be used in down stream analysis.

You may also notice that the order of the rows in adiExons$annotation is exactly the same as the order of the exons in the input GTF file.

Cheers,

Yang

Quick question, Im trying to run featureCounts with DEXSeq.

> head(adiExons$counts)
D104.bam D121.bam D153.bam D155.bam D161.bam D162.bam D167.bam D173.bam
WASH7P 0 0 0 0 0 0 0 0
WASH7P 9 4 2 6 4 5 5 7
WASH7P 8 1 2 5 0 1 0 2

How do I know which exon these are? Or more importantly, how does DEXSeq get the information it needs?

[sindrle]

Quick question, Im trying to run featureCounts with DEXSeq.

> head(adiExons$counts)
D104.bam D121.bam D153.bam D155.bam D161.bam D162.bam D167.bam D173.bam
WASH7P 0 0 0 0 0 0 0 0
WASH7P 9 4 2 6 4 5 5 7
WASH7P 8 1 2 5 0 1 0 2

How do I know which exon these are? Or more importantly, how does DEXSeq get the information it needs?

[sindrle]

Thank you for the clarification! Then I know what I did wrong, and maybe thats why it crashed..

[yangliao]

It is good that we have found the problem.

Every meta-feature in featureCounts can have one or multiple features. Each row in the GTF file is a feature if the third column equals the feature type you wanted (e.g., "exon", or any value defined in variable "GTF.featureType"). FeatureCounts then groups the features that have the same meta-feature-id into meta-features.

For example, if a row in the GTF file contains gene_id "LOC729737", then the meta-feature-id for this feature is LOC729737.

If the attribute name is not "gene_id" in the GTF file, you need to specify the correct attribute name. For example, if a row in the GTF file contains feature-name "LOC729737", you need to give a parameter telling featureCounts the new attribute name: GTF.attrType="feature-name", so that featureCounts can find the correct meta-feature-id by searching the attribute name (feature-name) in every row in the GTF file.

Thereby, if you want to assign reads on meta-feature level, please set GTF.attrType="gene_id", GTF.featureType="exon" and useMetaFeatures=TRUE.

If you want to assign reads on feature level, please set all variables the same values as above, but only useMetaFeatures=FALSE. FeatureCounts will do what it should do.

FeatureCounts counts reads for every gene in one single run. The returned R object contains a count table where the rows are meta-features (or features, depending on what you need) , and the columns are the input SAM/BAM file names. Every number in the table is the number of reads (or fragments if you specified isPairedEnd=TRUE) in that input file assigned to that gene.

Cheers,

Yang

I see!

R crashed making the "$counts" table.

featureGENE = "gene_id"

My goal was to count reads "at gene level", thus adding all exon reads for one gene together.

Is it, be the way, possible to count more than one feature in each run?

[sindrle]

The only problem is that your command line has an additional parameter "GTF.featureType=featureGENE". What was the value of variable featureGENE? The value should be "exon" for the GTF file you're using.

I see!

R crashed making the "$counts" table.

featureGENE = "gene_id"

My goal was to count reads "at gene level", thus adding all exon reads for one gene together.

Is it, be the way, possible to count more than one feature in each run?

[yangliao]

Dear Sindrle,

Sorry for the crash of R -- did it crash within a function of Rsubread or somewhere else? The ".featureCount" files are indeed very big -- they contain the assignment result of each single read (or fragment) in the input SAM or BAM files, and is uncompressed. That is why they are usually larger than the input BAM files.

The format of the GTF file seemed perfectly correct. I tried to put the first couple of lines you provided into a GTF file, and ran featureCounts on a paired-end BAM file. The result seemed correct -- featureCounts loaded 3 features and 2 meta-features from the GTF file, and assigned 286 reads to the 2 meta-features. The R object returned from featureCounts has the correct numbers on each meta-feature.

The only problem is that your command line has an additional parameter "GTF.featureType=featureGENE". What was the value of variable featureGENE? The value should be "exon" for the GTF file you're using.

Cheers,

Yang

        ==========     _____ _    _ ____  _____  ______          _____
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
       Rsubread 1.12.6

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           o A.bam                                          ||
||                                                                            ||
||             Output file : ./.Rsubread_featureCounts_pid24234               ||
||             Annotations : in.gtf (GTF)                                     ||
||                                                                            ||
||                 Threads : 1                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||                                                                            ||
\\===================== [url]http://subread.sourceforge.net/[/url] ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file in.gtf ...                                            ||
||    Number of features is 3                                                 ||
||    Number of meta-features is 2                                            ||
||    Number of chromosomes is 1                                              ||
||                                                                            ||
|| Process BAM file A.bam...                                                  ||
||    Assign reads to features...                                             ||
||    Total number of reads is : 1124980                                      ||
||    Number of successfully assigned reads is : 286 (0.0%)                   ||
||    Running time : 0.05 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
\\===================== [url]http://subread.sourceforge.net/[/url] ======================//

Im a bit disappointed, after finishing read summary after 5.5 hours, R just crashed.. R version 3.0.2 Rsubread version 1.12.6.

Maybe I can read the "reported output", the .featureCount files thats 2x the size of the BAMs (!)

Heres the GTF:

chr1 unknown CDS 69091 70005 . + 0 gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown exon 69091 70008 . + . gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown start_codon 69091 69093 . + . gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown stop_codon 70006 70008 . + . gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown exon 134773 139696 . - . gene_id "LOC729737"; gene_name "LOC729737"; transcript_id "NR_039983"; tss_id "TSS18541";
chr1 unknown exon 139790 139847 . - . gene_id "LOC729737"; gene_name "LOC729737"; transcript_id "NR_039983"; tss_id "TSS18541";

[sindrle]

Im a bit disappointed, after finishing read summary after 5.5 hours, R just crashed.. R version 3.0.2 Rsubread version 1.12.6.

Maybe I can read the "reported output", the .featureCount files thats 2x the size of the BAMs (!)

Heres the GTF:

chr1 unknown CDS 69091 70005 . + 0 gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown exon 69091 70008 . + . gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown start_codon 69091 69093 . + . gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown stop_codon 70006 70008 . + . gene_id "OR4F5"; gene_name "OR4F5"; p_id "P1230"; transcript_id "NM_001005484"; tss_id "TSS14428";
chr1 unknown exon 134773 139696 . - . gene_id "LOC729737"; gene_name "LOC729737"; transcript_id "NR_039983"; tss_id "TSS18541";
chr1 unknown exon 139790 139847 . - . gene_id "LOC729737"; gene_name "LOC729737"; transcript_id "NR_039983"; tss_id "TSS18541";

[yangliao]

Dear Sindrle,

Could you please provide a couple of lines in "UCSC hg19 genes.GTF"? There are some GTF or GFF files not in the format that is currently supported by featureCounts.

Yang

Hi!
Im trying your package, using UCSC hg19 genes.GTF.

Im getting this error:

|| Load annotation file genes.gtf ... ||
|| Number of features is 0 ||
|| WARNING no features were loaded in format GTF. ||
|| annotation format can be specified using '-F'.

This is my code:

featureCounts(files=BAMs,file., annot.ext=gtf,isGTFAnnotationFile=TRUE,useMetaFeatures=TRUE,GTF.featureType=featureGENE,GTF.attrType=attributeGENE,nthreads=8, reportReads=TRUE)

"BAMs" are a string character with the file names
"gtf" is the gene.gtf
"attributeGENE" is "gene_id"

EDIT:
Had to delete this: GTF.attrType=attributeGENE
guess it was unnecessary

[sindrle]

Here is my comparison so far:

HS 68,2 56,9 57,2 64,7 59,7 63,8 60,9
FC 62.5 56.8 57.1 64.6 59.5 63.6 60.7

HTseq was run with "intersection strict", so maybe thats why it counts a little bit more reads I guess.

I like that featureCounts is implemented in R, its easy to run
On my Macbook the running time is about the same as for HTseq, but maybe its because I had the option "reportReads = TRUE".

I also like that it reports gene length and may handle multiple input files, it also names the output according to the input automatically.

What will make or break it is how easily its implemented with the DEseq2 workflow.. HTSeq is very easy to use.

[sindrle]

Hi!
Im trying your package, using UCSC hg19 genes.GTF.

Im getting this error:

|| Load annotation file genes.gtf ... ||
|| Number of features is 0 ||
|| WARNING no features were loaded in format GTF. ||
|| annotation format can be specified using '-F'.

This is my code:

featureCounts(files=BAMs,file., annot.ext=gtf,isGTFAnnotationFile=TRUE,useMetaFeatures=TRUE,GTF.featureType=featureGENE,GTF.attrType=attributeGENE,nthreads=8, reportReads=TRUE)

"BAMs" are a string character with the file names
"gtf" is the gene.gtf
"attributeGENE" is "gene_id"

EDIT:
Had to delete this: GTF.attrType=attributeGENE
guess it was unnecessary

[shi]

Do not change your sam file, otherwise you may get unexpected results. Just change the 9th column of your gtf from for example

Name=exon1;created by=User;modified by=User

to

gene_id "exon1" (space delimited)

[chibouki]

I know, I stilled receive an error message about the column nine of the gtf (I may have changed it since I posted here for the first time...) but it worked online after cutting the sam