[shi]

I'm still not sure if it is the issue with paired-end reads that caused the problem. You can try to change those parameters to see it will work. But you may also try to count your reads as single-end reads by NOT using the '-p' option. This will tell us if the problem arose from dealing with the paired-end reads. Your command should be like this:

featureCounts -a mm10/annotation/mm10.allmrna.gtf -t exon -g gene_id -b -f -i tophat_out/accepted_hits.sort.bam -o subread_counts.txt

Wei

[choseqid]

Thanks, ddb, for the reminder. It looks like the problem lies in the way I aligned the reads with tophat: I allowed multiple hits (which apparently hampers the sorting by name) and didn't disable the separate alignment reporting for unpairable reads (ie. didn't use --no-mixed). Would fixing these two parameters help?

[ddb]

"featureCounts requires that for paired-end read data both ends must be included in the SAM/BAM file and the two reads from the same pair must be next to each other."

If this is not stated in the User Guide (I did not see it there) then it should be added as it is essential for correct functioning of the program.

[choseqid]

Dear Wei,

Thanks for the quick reply. Attached are the files you ask for. I am including the R output, as I do not have any from command line other than the error message I already quoted.

Cheers,
CHo

[shi]

Dear Cho,

Could you provide the complete output of your featureCounts run? It is hard to figure out what went wrong from the information you currently provided.

Also could you provide the first 100 lines of your annotation file and also the first 100 reads in your BAM file?

Cheers,
Wei

[choseqid]

Dear Wei,

I have been trying to run featureCounts, both from R and command line, but I keep getting a segmentation fault. I checked the different things suggested on this discussion thread and also tried to allocate more memory to the job resp. R session, but it didn't help. I use subread 1.3.6 and here is my command line:

featureCounts -p -P -d 50 -D 600 -a mm10/annotation/mm10.allmrna.gtf -t exon -g gene_id -b -f -i tophat_out/accepted_hits.sort.bam -o subread_counts.txt

Here is the error message:

/var/spool/gridengine/node-hp0211/job_scripts/1095023: line 10: 57569 Segmentation fault (core dumped)


Have I overseen anything?

Thanks in advance,
Cho

[shi]

Hi Bruce,

I agree it is better to count fragments for paired-end data. Looking at read counts is just to help diagnose what the problem was and it turned out this is quite helpful.

The percentage of multi-overlapping reads is much smaller that of multi-overlapping fragments, suggesting that something went wrong when read pairs were being summarized. We have never seen this before for the summarization of mapping results from Subread and a few other aligners.

One possibility is that in your mapping results the order of the two reads from the same pair was altered, ie the second read appeared before the first read. If this is the case, the read pair might be wrongly assigned. You may check a few multi-overlapping fragments to see if this is the case.

Alternatively, you may try other aligners as well. Subread is guaranteed to work with featureCounts.

Hope this helps.

Best wishes,

Wei

[bruce01]

Hi Wei,

I had looked at these before, but am keen to keep the counts to fragments as I believe there is more accuracy in this method. I have run for the initial trimmo sam, and the one with all non-pairs removed:

::::::::::::::
featco.pair.read.diags
::::::::::::::
35364121 ACCEPTED_GENE
4944952 MULTI_MAPPING
6848219 NOTFOUND_GENE
 338020 OVERLAPPED_GENES
::::::::::::::
featco.read.diags
::::::::::::::
35711285 ACCEPTED_GENE
4944952 MULTI_MAPPING
6902210 NOTFOUND_GENE
 340271 OVERLAPPED_GENES

[shi]

Hi Bruce,

I would suggest you trying to count reads instead of fragments to see if you will still get a large number of reads overlapping with multiple genes. This will help determine if the large number of multi-overlapping fragments you observed were due to summarization.

If you still get a large number, then that will mean it is either a mapping problem, or a problem with your data generation, or a lot of genes in your annotation overlapping with each other.

You can simply remove those paired-end parameters, such as -p -P -D and -C, from your command to summarize reads instead of fragments.

Best wishes,

Wei

[bruce01]

Hi Wei,

I was hopeful this was the issue resolved. Using Trimmomatic, you are returned 'paired' and 'unpaired' fastq files, so I immediately thought that was a nice way to remove unpaired reads; I only aligned the paired fastq. I have a script to remove all unpaired lines in SAM, and have confirmed the new, paired sam. When it is run in featureCounts, I still get the very large amount of OVERLAP:

::::::::::::::
trimmo.paired.diags
::::::::::::::
  87664 ACCEPTED_2VOTE_GENE
11663844 ACCEPTED_GENE
2824684 MULTI_MAPPING
1925844 NOTFOUND_GENE
7234680 OVERLAPPED_GENES
  10940 PAIR_DISTANCE
::::::::::::::
trimmo.diags
::::::::::::::
  87710 ACCEPTED_2VOTE_GENE
11713945 ACCEPTED_GENE
2826089 MULTI_MAPPING
1929703 NOTFOUND_GENE
7380972 OVERLAPPED_GENES
  10940 PAIR_DISTANCE

The difference is 201,703 counts, which is half the number of reads removed, ~70% are from the OVERLAPPED_GENES, but the rest are not 'correctly' allocated.

Any other ideas about this?

[shi]

Hi Bruce,

The three SAM files you sent to me contained different numbers reads, so I'm not sure if you used the same set of reads in your evaluations.

Also, your 'trimmo' SAM file is not in the correct format. One of the read pairs had only one read included in the file, the other was missing. The ID of this read pair is:

HWI-ST600:257:C1MUWACXX:3:1101:8762:2175

featureCounts requires that for paired-end read data both ends must be included in the SAM/BAM file and the two reads from the same pair must be next to each other. It has a unpredicted behavior if this requirement is not met. The problematic read pair in your 'trimmo' SAM file is located at line 17, so it is possible all the subsequent read pairs were incorrectly summarized. This may explain the bizarre result you observed for this dataset.

I speculate that the discrepancy in your previous results was also due to the missing reads in your SAM files.

Best wishes,

Wei

[bruce01]

Hi Wei,

so to return to my issue above with the many OVERLAP reads resulting from using Trimmomatic:

Having contacted the authors of STAR we have determined that it is not the aligner at fault for the issue I have reported to you. I have attached 3 short SAM files; I took the first 100 OVERLAP reads from Trimmomatic trimmed reads as determined by featureCounts. I grepped these out of Trimmomatic qulity trimmed at Q30 ('trimmo'), hard-trimmed by 10bp ('trim') and not trimmed ('notrim') SAM. I then used featureCounts to again determine what features these reads map to. Reads were initially found using:

grep -m 100 OVER trimmo.featco.counts.reads | cut -f 1 | grep - ./trimmo/Aligned.out.sam > trimmo.100-over.sam

now when I rerun featureCounts on the 'trimmo' set:

~/bin/subread-1.3.5/bin/featureCounts -a Bos_taurus.UMD3.1.70.gtf -i trimmo.100-over.sam -o trimmo.100-over.counts -p -P -C -R -D 500000

then

grep -v OVER trimmo.100-over.counts.reads | wc -l

returns 10. So now, from exactly the same 100 reads that were initially OVERLAP, 10% are assigned with clarity to a single feature. The other two sets, 'notrim' and 'trim', return the same results as previously, each having only 2 reads return OVERLAP.

Interestingly, or unfortunately, depending on viewpoint:

grep -f trimmo-100.over.sam notrim-100.over.sam | wc -l

returns 71. So 71 of ~200 lines are identical between 'trimmo' and 'notrim', ie there has been no trimming done on these reads and they align to the same position. But when these reads are under the header of 'trimmo' they produce OVERLAP, and under the header 'notrim' they produce ACCEPTED.

I have now taken all reads which are called OVERLAP in the trimmo analysis and grepped them from trimmo and notrim SAM, then compared those two files for number of lines which are identical.

wc -l trimmo.featco.reads.solo
7380972 trimmo.featco.reads.solo
wc -l notrim.trimmo-overlap.matched.lines.sam
5036198 notrim.trimmo-overlap.matched.lines.sam

So that's 68% identical, of which 100% in trimmo data are called OVERLAP, but those 5,036,198 from notrim data, identical reads to those called OVERLAP, give a breakdown thus (obviously there are nearly double the number of reads as this is PE data and the featureCounts run include -p -P flags):

28587 ACCEPTED_2VOTE_GENE
2616603 ACCEPTED_GENE
6191 NOTFOUND_GENE
55188 OVERLAPPED_GENES
234 PAIR_DISTANCE

Which is 96% ACCEPTED, 2% OVERLAP versus 100% OVERLAP from the trimmo data. Which is a massive 'discrepancy'.

Your thoughts on this would be greatly appreciated.

[Nicolas Nalpas]

Dear Wei,

That is great, very helpful.
Thanks a lot.
Best wishes,
Nicolas

[shi]

Dear Nicolas,

No worries. Thanks for your interest in featureCounts.

featureCounts calls an overlap when there is at least 1bp overlap found between a read and a feature. This also applies for paired-end reads, ie. an overlap is called as long as 1bp overlap was found between one of the two reads and the feature.

featureCounts will check every feature to see if it overlaps with a read. If more than one overlapping feature was found, it will not assign this read to any feature (by default). However if '-O' is specified, featureCounts will assign it to all its overlapping features.

For paired-end reads, featureCounts also checks if both ends from the same fragment overlapping with the features. For example, if a fragment has only one end overlapped with a feature F1 but this fragment also has both ends overlapped with a feature F2, this fragment will be assigned to F2 because this overlap has a stronger support.

Hope this makes it clearer. Don't hesitate to ask me if you have any further questions.


Best wishes,

Wei

[Nicolas Nalpas]

Dear Wei,

Really sorry to bother you again, I have just one question as a matter of interest regarding featureCounts.

How do you handle your read count summarisation, when is a read assigned to a feature exactly? I know that HTseq-count has three summarisation modes (see pictures enclosed): union, intersection_strict, intersection_nonempty; so I am just interested to know what is yours (sorry I was not able to find the information in the users guide and paper).



Thanks a lot for all your help.
Best wishes,
Nicolas