This is a test version of Biostars. For the public version, visit https://www.biostars.org.
featureCounts cDNA vs DNA files

featureCounts from the subRead package works quite well with cDNA reads aligned to genomic DNA_ref.fa files, but I am having difficulty getting featureCounts to work with cDNA files aligned to the cDNA_ref.fa.

mRatBN7 GTF File:

DNA ref.fa File:

cDNA ref File:

minimap2 used for mapping Nanopore cDNA reads over to the different reference files. The data revealed by samtools flagstat indicates that my cDNA data aligns better to the cDNA?_ref.fa file, and would like to use this data... if only featureCounts would cooperate. But I always get 0.00% reads aligned with the cDNA .

featureCounts( DNA.sam / cDNA.sam, annot.ext = GTF_mRatBN7, isGTFAnnotationFile = TRUE)

It was suggested to me to use samtools idxstats and using the table directly from there as my gene output, but it doesn't quite feel right.

Any suggestions?

PS. I have to use R for subReads package as it is not yet fluid with the newer Apple M1 chips, however I do not suspect this is the problem as I experienced this also on Ubuntu terminal.

cdna featurecounts subread rna
• 3,393 views • link
written by camppatrick • camppatrick

featureCounts can also be run as a separate tool, but as GenoMax stated, even the R version should run seamlessly on M1 chips

• 0 views • link
written by Istvan Albert Istvan Albert

1 answer

Problem in this case is likely that your GTF file references the chromosome names where as your cDNA file has names for cDNA's. So those two don't match. You could create a simple SAF format file for your cDNA simply listing the names of cDNA's and their start/stops. You could do that using the suggestion you got

It was suggested to me to use samtools idxstats and using the table directly from there as my gene output, but it doesn't quite feel right.

SAF file format example :

GeneID      Chr Start   End Strand
497097      chr1    3204563 3207049 -
497097      chr1    3411783 3411982 -
497097      chr1    3660633 3661579 -
100503874   chr1    3637390 3640590 -
100503874   chr1    3648928 3648985 -
100038431   chr1    3670236 3671869 -

I have to use R for subReads package as it is not yet fluid with the newer Apple M1 chips

Not sure what this means. Anything coded for intel chips should work on M1 via rosetta translation.

• 0 views • link
written by GenoMax GenoMax

Header format from bam file

Created SAF file

Hello, Thanks for the input. I created a SAF file of this format using the GTF file and some regex magic in R. I now am running:

featureCounts(<data.bam>, annot.ext = <saf_file.txt> )

And featureCounts runs immediately, returning 0.00% reads.

Things I have tried with no success:

  • saving the saf_file.txt as saf_file.saf
  • Putting SN: in front of each transcript name for the saf file
  • leaving "chr" out of the Chr column
• 0 views • link
written by camppatrick • camppatrick

Subread

Attached error message from trying to make the Subreader file for MacOS. Rosetta is installed on my computer, checked by forcing some apps to open using rosetta.

But I believe this issue is probably better suited for its own thread.

• 0 views • link
written by camppatrick • camppatrick

I created a SAF file of this format using the GTF file and some regex magic in R.

This is not correct. cDNA entries in your transcript files are standalone sequences. They are of length indicated in first screenshot by LN field. You can't substitute chromosomal coordinates in your SAF file. So the entry in your SAF files should look like this.

ENSRNOT00000033201.6 chr12 1 276 -

When you count the BAM file aligned against the cDNA reference these are the sequences/lengths you need to use.

Edit: Added chr field that was missing.

• 0 views • link
written by GenoMax GenoMax

The subread user guide says these columns are necessary : GeneID Chr Start End Strand

Leaving out the Chr column generates the error that the data is not in saf format.

Unfortunately, adding in the chromosome column (and playing around with the formatting) still results in 0.0% alignments. Am I just on the wrong path here? Should Nanopore cDNA reads just be mapped against the genome and not the transcriptome? Programs such as IGV, htseq, and featureCounts all seem to work a lot smoother with the genomic mapping, and there really aren't any forum threads online of people facing this problem.

• 0 views • link
written by camppatrick • camppatrick

My apologies for missing that.. Can you add chr names from your other file but use the lengths of transcripts as start and stop?

• 0 views • link
written by GenoMax GenoMax

saf file

featureCounts output

I am posting the example saf file I made. To note is that I tried many different variations of formatting on the GeneID, Chr, and Strand columns, ie: changing +/-, with/without "chr", etc... I also made sure to include transcripts that had high counts in the idxstats output.

These all yield 0% reads.

• 0 views • link
written by camppatrick • camppatrick

Keep in mind that featureCounts does not count multi-mapping reads by default. If you have more than one transcript for genes, this is likely happening. You can add -M option to count multi-mapped reads.

In any case you transcripts are small and not sure how long your nanopore reads are. Examining alignments carefully is in order.

• 0 views • link
written by GenoMax GenoMax

Hm. Unfortunately I don't think it is an issue of multi mapping, but still rather featureCounts is not interpreting some syntax somewhere.

The option -M also produces 0 reads.

To note: I am using a small subset of data for quick alignment and pipeline testing. Hence the low number of reads (this is also just barcode 1 of 12. My reads are 100-1000bp in length, and I have a few thousand reads mapping to the genome (minimap -ax splice ref = genomic file), and ~30000 total primary alignments to the cDNA reference.

Perhaps another column of information is needed in the .saf file? Or a specific syntax of some column is false. I am at a loss.

• 0 views • link
written by camppatrick • camppatrick

Log in to answer this question.

More posts like this