I have a set of scRNA-seq samples enriched with FACS for cells expressing a specific gene reporter (TdTomato). In particular the gene I want to report has positive counts in the resulting matrix for 97% of cells.
I followed CellRanger documentation to create a custom reference. I added TdTomato on an additional chromosome and I added an entry in the gtf:
TdTomato . exon 1 1431 . - . gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato"
The command I'm using is:
cellranger-2.1.1/cellranger mkref --genome=refdata-cellranger-mm10-2.1.0.tdtomato --fasta=genome.tdtomato.fa --genes=genes.tdtomato.gtf
genome.tdtomato.fa is the original reference genome (mm10, as provided by cellranger) with an additional chromosome at the end with header:
genes.tdtomato.gtf the original reference genome gtf with the additional line indicated previously at the end of the file.
I actually tried both on the '+' and on the '-' strand, however I guess this should not affect the read mapping, since TdTomato is on an independent chromosome.
cellranger count using the custom reference:
s="SAMPLE" cellranger-2.1.1/cellranger count --id="$s"_tt --transcriptome=refdata-cellranger-mm10-2.1.0.tdtomato --fastqs=original_fastqs_directory/"$s"/fastqs --sample="$s"
I am using the default parameters of CellRanger.
However, I found that less than 1% of cells seems to express TdTomato. If I use the '+' strand I get slightly more counts (2.4% of cells have positive TdTomato).
What I checked:
- With IGV I can see that there are reads mapping to TdTomato
- Using samtools I am able to see the same number of reads mapped on TdTom either with the + or - strand
- However, following this documentation I noticed that by executing the CellRanger filtering steps, more reads are filtered in the case of the - strand
This reduces the number of reads used for counts to 44 (2.4% of the original 1807 reads, - strand), 240 (13.3%, + strand) for TdTomato, 71620 (21.2%) for the gene reported by TdTomato. Interestingly, there are fewer reads for TdTom than counts, so there may be some kind of normalization afterwards?
The gene reported by TdTomato is longer, 11394 bp than TdTom, 1431 bp. The total length of the exons for the two genes is 5349 vs TdTom: 1431. Does this influence exponentially the amount of reads mapped to the two genes?