samtools to count the number of reads mapped to each spike-in for each sample

Question

[kstangline]

My goal is to use STAR to create a new genome with the spike-ins listed below by combining both hg38.fa and spike-in. Once I have the genomes created, I'll align FASTQs to this newly created genome and use samtools (or another tool) to count the number of reads mapped to each spike-in for each sample. Would anyone know how to use samtools in this case?

Spike-in 150bp:

AGGCTTTGTTATTGGTATTGACTTACAAACAGTTAAGCCATTTGAATATGATAATGTAGTTGCAATAAAAGGAGATTTCACCTTAGAAGAAAATTTGAACAAAATTAGAGAGCTAATTCCAAATGATGAAAAAAAGGTGGATGTGGTTAT

Spike-in 250 bp:

ACGATAAAATTGGTTTTGCCTTTCAGCAATTCAACTTAATTCCTTTATTAACTGCCTTAGAAAATGTTGAACTTCCACTGATTTTTAAATATAGGGGAGCAATGAGCGGAGAAGAGAGGAGGAAGAGAGCTTTAGAATGCTTAAAGATGGCAGAGTTGGAGGAGAGATTTGCCAATCACAAACCAAATCAGTTGAGTGGAGGGCAACAACAGAGAGTTGCTATAGCGAGGGCTTTGGCAAACAACCCACC

Spike-in 350 bp:

CACCTGTCACATTTCCAATCGGCTCCAGGAAGAGAGAAGTGACGGCTTGATCCTGTAGTAATCCGGGATCGACTTAAGGGGTGCAGCGACCACGGCGGATCGGGCGTCGCAATAGTCCTCCTGTTAGGAGGGTCCTTCTAATGTTAACGCCCGAATATTAGTCATATTTTGCTAGCGCCTATCAGCGTAAGATATGATTTAAGTTACACCAGGAGAGTAGCGAGATAGAACCACTCGTTGGATCGGTCTTTCTTAATTGACTACTATCAGATCCGGCGCATGGCGCTGAGGTCAAACTACATTACAGGCCCTGGTTTCCATGGGTCAGCGCAAGTACAGGCGAGCAGATA

Spike-in 450 bp:

CACAAGAATCCCTGCTAGCTGAAGGAGGGTCAAACTATAACACCTTTAGCATTCGTACAGGCAGGCTAAGTGAATACTAACCCACCGGCAGCCCGTTGTAGTAACGTTGACCCCTGGCTCGGAGACATTTGGTGTTGCCTAGTACTAGGTGACTGGTACCGATTCATAGGTTCGCCATTCTCTTATCGAGAGCCCGAGGTAGACTATCTTCCAGATGATGCCATACGTTCACTCAATCGCGCGGCATGCACGGTGGGGCTACGAACTTGCTATCCATAGGCTCTAGATGTGGTAGAAATATGCTGCAGGGGTTCTGTCGAATTTGCTCGGCAACCGTGGCCGTGTATGCTTTCATATCCCGGCGGTGTGATCTAGCCTTCTCGCCATATGAGGGCGCTGAGCATAGACCCAAACCCGACTAGTCGAATCTTAGGGTTGTATGCTAGAACG

[GenoMax]

Add the spike-in's as individual references in your genome file. After alignments, use samtools idxstats to get reads aligned to each of the spike-ins from BAM.

Answer 1

[Istvan Albert]

You can do what GenoMax suggested but I would recommend an alternative that fits with RNA-Seq workflows. Use an annotation information usually in the form of the GFF file that comes with your spike data (or create a file like that) .

Now treat your spike as if it were another gene on the genome and incorporate it into the entire process. You'll be able to validate the statistical process's powers that way.

Basically instead of samtools use a tool like featureCounts to count the reads that span over a region. The counting tool is designed to count reads and will provide you with additional utility that you will likely need later on.