[GenoMax]

That's really strange. It is, indeed, not a valid sam file.
All of the columns are present. I have no idea how you were able to generate output missing some columns... could it be some kind of copy-paste error? Or did you reformat the data in some way?

When I copied and pasted the data last night something must have gone wrong. I have recopied the data in post above. Sorry about that.

[Brian Bushnell]

That's really strange. It is, indeed, not a valid sam file. I just did some mapping and generated this:

@HD	VN:1.4	SO:unsorted
@SQ	SN:ecoli_K12	LN:4639675
@PG	ID:BBMap	PN:BBMap	VN:35.51	CL:java -ea -Xmx1g align2.BBMap in=80x_perfect.fa.gz reads=1000 int=f out=xx2.sam ow
400	16	ecoli_K12	778201	45	150=	*	0	0	TCTTCCGGCAACTGATGGACAGGTCAAATTCCCTGCCTGGTCGCCGTATCTGTGATAATAATTAATTGAATAGTAAAGGAATCATTGAAATGCAACTGAACAAAGTGCTGAAAGGGCTGATGATTGCTCTGCCTGTTATGGCAATTGCGG	*	NM:i:0	AM:i:45
401	0	ecoli_K12	940032	45	150=	*	0	0	GTAATACTACTTTCGAGTGAAAATCTACCTATCTCTTTGATTTTCAAATTATTCGATGTATACAAGCCTATATAGCGAACTGCTATAGAAATAATTACACAATACGGTTTGTTACTGGAATCAATCGTGAGCAAGCTTGAGTGAGCCATT	*	NM:i:0	AM:i:45
402	0	ecoli_K12	543578	45	150=	*	0	0	CATTTCATCGCAAATCGCCCGCATGTCGTTTTCTAACTGTTGGGTGAAATCGCGCAGCACGGCAGCGTCGGTATGACGACAATCAATGGTGAACGTGGTTTTACCCGGCACCACATTTACCGTATTCGGGCGCGGCTCTACTTTGCCAAA	*	NM:i:0	AM:i:45

All of the columns are present. I have no idea how you were able to generate output missing some columns... could it be some kind of copy-paste error? Or did you reformat the data in some way?

[shi]

Hi GenoMax,

Thanks for providing the sample reads. However, we noticed that there are missing columns in your mapping results (either PNEXT or TLEN or both). This might cause problems for read counting by featureCounts and lead to incorrect counts. For example, featureCounts might not be able to find NH tags which are used to identify multi-mapping reads.

[GenoMax]

Hi Shi,

Brian is on the right track.

Here is a sample of reads

With v.1.3 tags

XXXXXXX:XXX:XXXXXX:2:1213:18683:18437 1:N:0:        0       chr10   376     3       50M     *       0       0       GTCTGTGTTCTGTGTCAGAGGAGCAGTAGAGCTGAGGCTATCAAGCTATC    BB@BBHEHHCG@<FC@CE@FGH@EGHGFEEHIHFEFE?GEHHHGEHHFHE       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1203:10888:83435 1:N:0:        0       chr10   413     3       50M     *       0       0       CTATCAAGCTATCAGGACATCCCAGATGCCCAAGACTACACTAGTCACCA    DDDBD?HFEE?GHGHEEEFIFIIHGHHHI?HIIIEHIIIIEHHHEHFHHI       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1114:19668:64532 1:N:0:        0       chr10   938     2       5M3I42M *       0       0       CCACTGACACCATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGT    DDDDDIIIIIGIIIIIIIIIHIIIIIIIIHIIIHIIIIIIIIIIIHIIIF       XT:A:R  NM:i:4  AM:i:2
XXXXXXX:XXX:XXXXXX:2:1103:21004:71100 1:N:0:        0       chr10   946     3       50M     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    0<D0@@FHHIIIIIE?F@GHE?H@GHHIIIEGH?HFE@GHHHH?ECHHII       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:2106:11030:47701 1:N:0:        0       chr10   946     3       50M     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    DDDDDHIFIIIHIIIIIIIIHHEHEHIIGHHIIIIHIIHHIIIHHHHHHE       XT:A:R  NM:i:0  AM:i:3

With v.1.4 tags

XXXXXXX:XXX:XXXXXX:2:1213:18683:18437 1:N:0:        0       chr10   376     3       50=     *       0       0       GTCTGTGTTCTGTGTCAGAGGAGCAGTAGAGCTGAGGCTATCAAGCTATC    BB@BBHEHHCG@<FC@CE@FGH@EGHGFEEHIHFEFE?GEHHHGEHHFHE       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1203:10888:83435 1:N:0:        0       chr10   413     3       50=     *       0       0       CTATCAAGCTATCAGGACATCCCAGATGCCCAAGACTACACTAGTCACCA    DDDBD?HFEE?GHGHEEEFIFIIHGHHHI?HIIIEHIIIIEHHHEHFHHI       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1114:19668:64532 1:N:0:        0       chr10   938     2       2=1X2=3I42=     *       0       0       CCACTGACACCATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTC
AGT     DDDDDIIIIIGIIIIIIIIIHIIIIIIIIHIIIHIIIIIIIIIIIHIIIF      XT:A:R  NM:i:4  AM:i:2
XXXXXXX:XXX:XXXXXX:2:1103:21004:71100 1:N:0:        0       chr10   946     3       50=     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    0<D0@@FHHIIIIIE?F@GHE?H@GHHIIIEGH?HFE@GHHHH?ECHHII       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:2106:11030:47701 1:N:0:        0       chr10   946     3       50=     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    DDDDDHIFIIIHIIIIIIIIHHEHEHIIGHHIIIIHIIHHIIIHHHHHHE       XT:A:R  NM:i:0  AM:i:3

Edit: Just saw the last post from Wei. Will look forward to a future release of subread.

Edit2: There was an error in the copy/paste of the data in original post. I have replaced the original with correct copy.

[shi]

featureCounts currently does not support letters '=' and 'X' in the CIGAR string and we think that is the reason why it failed to assign most of the reads in your v1.4 data. We will add support to these letters in featureCounts.

We are going to make a major release of our software package next week but unfortunately we are not able to include this support in the release. We will work on fixing this issue after the release. Will let you know once this is done.

[Brian Bushnell]

For what it's worth, the only difference between SAM 1.3 and 1.4 format output BBMap is that in 1.3 mode, it use the "M" symbol in the cigar string; in 1.4 mode, it uses "=" and "X" for match and mismatch in the cigar string. 1.4 will still occasionally contain "M" symbols when there is a "N"-called bases, as that does not strictly speaking match or mismatch. Given that a small fraction of the reads do get assigned from v1.4 sam files, I wonder if featureCounts might be using only the reads that contain at least one M operation in the cigar string?

[shi]

Thanks @GenoMax,

Could you also provide mapping results for a few reads for each version?

[GenoMax]

Hi @GenoMax,

featureCounts does support SAM specification version 1.4. Have you looked at the .summary file generated by featureCounts to see why reads were not assigned?

Best,
Wei

Hi Wei,

Here is a comparison of the two flags for a sample. I have used the -M flag for featureCounts.

With SAM v.1.3 flags

Assigned	     104867565
Unassigned_Ambiguity	577037
Unassigned_MultiMapping	0
Unassigned_NoFeatures	52436519
Unassigned_Unmapped	1772756

With SAM v.1.4 flags

Assigned	      66620
Unassigned_Ambiguity	280
Unassigned_MultiMapping	0
Unassigned_NoFeatures	157814221
Unassigned_Unmapped	1772756

[shi]

Hi @GenoMax,

featureCounts does support SAM specification version 1.4. Have you looked at the .summary file generated by featureCounts to see why reads were not assigned?

Best,
Wei

[GenoMax]

Hi Wei,

Do you have plans to add support for SAM v.1.4 format tags to featureCounts?

We use BBMap for alignments and it outputs SAM v. 1.4 tags by default. We have been having trouble assigning read counts using those alignments with featureCount. featureCount does not generate an error but only assigns a small number of (< 1%) reads.

Note: BBMap can output SAM v.1.3 tags and then the counting works, which we are using as a workaround for now.

Thanks.

[shi]

Hi Bruce,

I agree with Devon that if your library is stranded you should count reads at the stranded mode (s=1 or s=2) rather than unstranded mode (s=0).

Also as I said earlier, my concern with performing stranded counting for unstranded reads is that you may have reduced capability to detect reads that overlap with more than one gene. This may result in reads being assigned to wrong genes.

Another thing you need to consider is that the annotated strands for genes might be incorrect, particularly for those computationally predicted genes. Stranded read counting will be problematic for such genes.

Wei

[bruce01]

The incorrect one will have vastly lower numbers.

This is my problem, I have samples where s=1, which is incorrect based on the library prep protocol which was dUTP based, have more counts than s=2. This is from FFPE tumor samples, and so I expect it to be bad, but not like this. Will go with s=2 anyway, thanks for your help.

[dpryan]

If your library is stranded, then use either s=1 or s=2, as appropriate. It seems that s=2 would typically be correct, since most stranded kits use a dUTP-based method (so the sense strand is the opposite of how read #1 aligns). When in doubt, featureCounts is really fast, so you can run it with s=1 and s=2 separately and then just determine which one was the correct setting from the numbers output (afterall, sometimes you don't know the strandedness). The incorrect one will have vastly lower numbers.

[bruce01]

Strand refers to sense and antisense, not +/-.

Sorry about the nomenclature, I have +/- in my head as they are in GTF and that is how I denote it in there.

Can I ask your opinion on what way to run featureCounts with stranded libraries? Presumably if strand is known, then using both s=1 and s=2 gives the more accurate counts that using s=0, particularly for overlapping genes, when strand is known?

[dpryan]

Strand refers to sense and antisense, not +/-.