Could someone explain why the default pipeline uses both STAR for alignment and then use the BAM files through Salmon for generating counts tables? Why not use featurecounts?

This was ported from the #CE8SSJV3N workflow I believe. Salmon handles transcript level quantification better than featurecounts. That being said, I am not sure how I feel about applying salmons probablistic assignment of reads to Ribo-Seq data. The coverage that Salmon expects in RNA-Seq will not be seen in Ribo-Seq. This probably warrants investigation

featurecounts for rnaseq data is very old-school. multi-mapping reads are a huge problem (they're just chucked out) and the presence of gene families throws things off in ways that are impactful for e.g. human disease. I think some of that reasoning may apply quite well to riboseq data. That said, it is a very special case given the short reads, so there are probably some optimisations to be made.

Interesting...so are these alignment stats not reflective of the salmon counts table? I assumed STAR aligns all the RPFs, throwing out any "Secondary Alignments", thus for SRX117080885 the initial 77.9M reads gets chopped down to only 13.1M reads, 16% of the original read count to be aligned and used for all downstream processes

Ribo-seq analysis isn't standardised because the questions people look to answer using Ribo-Seq are not standardised. Depending what you wish to investigate you will use different steps each with their own sets of pros/cons. Pretty sure Salmon will take reads aligned to the transcriptome so if these are genome alignments then you would expect some differences yes. Apologies, I'm not quite clear on your question so I doubt I have answered it

I appreciate your insight and time to answer my questions--to clarify: From the pipeline it seems that STAR assigns the majority of the reads as Secondary Alignments (which I understand as multi mapped reads). I believe these are not included in the BAM file that gets fed into Salmon, and only the 13.1M reads that were uniquely mapped. Is this a correct understanding? I ask because @Jonathan Manning said that Salmon > than featureCounts because it does use the multi-mapped reads (vs just chucking them out). So does the 13.1M aligned reads in the Alignment Panel of MultiQC not accurately reflect the amount of reads being used for Salmon.merged.counts files? And along those lines, if only 16% of total reads are useful after aligning to an index created from the whole primary assembly, would it be incorrect to generate an index from the Homo_sapiens.GRC38.cdna.all.fasta? As my colleague says they retain 50%+ of their total reads.

Ambiguous mappings can also mean that the one alignment could be assigned to two different features. These are the ones that featurecounts throws away but Salmon does not. The BAM given to Salmon should have secondary alignments but Salmon does not consider them as far as I know.

Thanks for that clarification. If I understand correctly, STAR Aligner will first generate a BAM file with 13.1M reads, where these reads may align at multiple locations within one gene. The other 64.8M reads that were assigned as "Secondary Alignments" were reads that were similar to the primary alignment in nucleotide sequence, but were mapping to many other genes. These are not passed onto Salmon. However, the 13.1M reads are passed onto Salmon, and a proportion of these reads may multi-map the same gene at different loci (e.g., exon 1, exon 5, exon 10), and Salmon uses probabilistic assignment to choose just one loci?

Dont mix up reads and alignments. One read can have up to 20 alignments if you used the default args. You give trimmed reads to STAR. It aligns each read returning alignments in BAM format. This BAM also has the secondary alignments so they are given to Salmon. But Salmon (to the best of my knowledge) will see if each alignment is a primary or secondary and only count the primaries. The example with exon1, 5 and 10 would all be secondary alignments as these come from different genomic loci. What Salmon does is when 10 transcript isoforms contain exon 1 it probabilistically distributes its "weight" across the possible transcripts using expectation maximisation.

Okay, but if that's the case then why does samtools stats output 13.1M total sequences, vs samtools flagstat that stated there were 77.9M reads total, with 64.8M secondary alignments (and the difference being 13.1M)? The way I interpret that is that 13.1M total sequences/reads are being passed through STAR, with 64.8M secondary alignments being filtered out.

samtools stats

ignores secondaries and

samtools flagstats

deals with all alignments but breaks them down https://www.biostars.org/p/412621/

https://nf-co.re/riboseq/1.0.1/results/riboseq/results-3d3e57978c5e0a418e56728496b6cc12700e6443/quantification/salmon/SRX11780885/logs/ Interestingly, the salmon_quant.log file states even less reads...?

hey, if i may chip in... I generally find it way easier to check the STAR section of the multiQC directly to interpret read fate at this stage. So from this you can see that only ~17 ish million reads went as input to STAR for SRX117080885. The data filtering flow is as follows (by default): input reads -> filter out too short after trimming adapters and low quality bases (here you are left with 46 mil reads) -> filter out rRNA contamintants etc with sortmeRNA (here over 20 mil reads map to contaminants) -> map remaining reads - ~17ish mil with STAR, where anything unmapped will also be filtered out. so you're left with ~13mil mapped reads. By default STAR maps to the genome but will also output a BAM with transcriptomic coordinates. if any UMIs, then reads will also be deduplicated, so you might get even less useful reads. The dedup transcriptomic BAM is used with salmon, and multimapped reads are considered.

Thanks for your insight--diving deeper into this I find that I have more questions...is the input reads number for -885 77.9 million? The fastqc for raw reads should reflect this, but it displays ~47 million. After these reads are trimmed, it shows ~46 million sequences/reads. Am I missing something?

no you're right, input is 47 million based on raw fastqc, i dont know where the 77.9 comes from, can't trace this right now

Perhaps the rawfastqc is bugged? It might be 77.9million to start, as Cutadapt Filtered Reads shows the same number as the FastQC (raw)

I guess it is because of the multiqc column label but flagstat does not handle reads. It handles alignments. Thats 77.9 million alignments

77.9 million alignments, meaning 47 million reads map to 77.9 individual positions?

yes @Jack Tierney, that's prolly it! the 77.9 mil is also present in the general statistics the very first section of multiQC so I can see how one can find that confusing..it's labelled "reads mapped"

@Ira Iosub you will know better than me, is this something we should change in #C05785S3TMX MultiQC usage or is it an issue with MultiQC?

need to think about it a little, and maybe check some multiQCs for rnaseq with samples I know well

@Richard Sam Yes and some of those dont map anywhere

might be sth general with just how the general stats table is built, so I will check in some rnaseq reports i have and let you know

Should the raw FastQC, to the Cutadapt reads passing filters, to the trimmed FastQC all be <2%? And somewhat unrelated--if I don't have UMIs, is it conventional to use SAMtools dedup to filter out PCR duplicates?

Are you saying only 2% of the reads pass cutadapt? 2% of the reads in the raw FASTQC make it to the trimmed FASTQC? How would you know which reads are duplicates if you don’t have UMIs?

My understanding is that SAMTools marks reads with the same starting sequence and end sequence as a duplicate. While it seems more useful for RNAseq, I wasn’t sure if it’s conventional for Riboseq given identical RPFs could be biological? https://www.htslib.org/algorithms/duplicate.html

I definitely wouldn’t advise it. It would be a detrimental to drop all of those kinds of duplicates

It seems for the test dataset, the # of reads in FastQC raw are very close to exactly the same number in Cutadapt, and just a slight decrease in the FastQC trimmed.

hey @Jack Tierney I looked at some multiQC reports of RNA-seq with

rnaseq 3.18

, but the multiQC looks quite different from this one in the general stats section ... and they seem to be reporting read counts there rather than alignments (not 100% sure). Specifically for riboseq we need to check the logs and where the multiQC is getting the info from. Slightly unrelated, even in the rnaseq multiQC reports there are a lot of columns with the same name "seqs" and unclear what they actually represent without chasing samples in the report. e.g. had to infer the original number of reads in the raw FASTQ by navigating to the fastqc (raw) section. or there are multiple "M aligned" columns - it's perhaps worth looking into making these tables more intuitive or clearly labelled for both rnaseq and riboseq, and potentially add checks that the numbers after various seps do add up to the previous step/original