SMART-seq2 单细胞转录组分析workflow

0. 质控

这一步主要看接头，其他由于 ERCC 等等原因，会和 fastqc 提供的正常报告有区别。

0.1 使用 fastqc 批量质控

find ./RAW/ERR*fastq.gz | xargs ./app/fastqc -t 6 -o ./fastqc_result

0.2 使用 multiqc 整合 fastqc 的结果

cd fastqc_results
multiqc .
firefox multiqc_report.html

• 用firefox之前的登录记得要用 ssh -X

1. 比对

如果有 spike-in，参考基因组应该加上 spike-in 的序列，同时，gtf 文件也要加上 spike-in。
spike-in fasta和gtf的下载地址

cat ERCC92.fa >> GRCh38.fa
cat ERCC92.gtf >> GRCh38.annotation.gtf

1.1 比对并查看比对质量

nohup ls *gz | while read id; do STAR --runMode alignReads --runThreadN 4 --readFilesIn ${id} --outFileNamePrefix ../alignment_result/${id%%.*} --genomeDir ../ref --readFilesCommand zcat; done > alignment.log 2>&1 &

# 整合多个文件的总reads数，比对率等等
ls *.final.out | while read id; do (cat ${id} |sed '1i ID |'${id%%.*} | awk 'BEGIN{FS="|"} {i=2;while(i <= NF){col[i]=col[i] $i " ";i=i+1}} END {i=1;while(i<=NF){print col[i];i=i+1}}' | sed 's/[ \t]*$//g' >> QC.matrix);done

1.2 生成bam文件并索引

samtools view -@ 10 -S SRR0011223344.sam -b > SRR0011223344.bam
#.bam排序，默认按照染色体位置
samtools sort SRR0011223344.bam -o SRR0011223344.sorted.bam
#索引
samtools index SRR0011223344.sorted.bam

2. 统计 countmatrix

2.1 使用 featrueCounts

../app/subread-2.0.0-Linux-x86_64/bin/featureCounts -T 6 -t exon -g gene_id -a ../ref/GRCh38_ERCC.gtf -o CountsRaw.txt  *.sorted.bam 1>FeatureCounts.log 2>&1 

3. 使用 scater 进行下游分析

不同的 r 对应安装的 scater 版本不一样。最新版本的 scater，需要 r 4.0。且不同版本的 scater 用法有一点不同。这里用的是 scater 1.15.32 。

#当前版本的测试数据
example_sce = mockSCE()

3.1 导入数据并生成 SingleCellExpriment对象

#FeatureCounts得到的matrix前五列都是基因信息，第六列之后才是counts。
countmatrix2=countmatrix[,-(1:5)]

#把列名改一下，FeatureCounts的输出是以文件名为列名，带有 sorted.bam 后缀
y=sapply(colnames(countmatrix2),function(x) strsplit(x,split = '.sorted.bam'))
colnames(countmatrix2)=y

#读取细胞信息，行名为细胞，列名为细胞的各种信息
cellinfo=read.table("./cellinfo.txt",header=T,sep = '\t',row.names = 1)

#as.matrix这一步很重要，不然后续会报错类型不对
sce=SingleCellExperiment(assays=list(counts=as.matrix(countmatrix2)),colData=cellinfo)

3.2 标记 spike-in 和 线粒体基因

is.spike <- grepl("^ERCC-", rownames(sce))
sce <- splitAltExps(sce, ifelse(is.spike, "ERCC", "gene"))

is.mito=grepl("^MT", countmatrix[1:nrow(countmatrix),1])

3.3 进行自动质控

sce<- addPerCellQC(sce, subsets=list(Mito=is.mito))
reasons <- quickPerCellQC(sce, percent_subsets=c("subsets_Mito_percent",
    "altexps_ERCC_percent"))
colSums(as.matrix(reasons))
sce$discard=reasons$discard

3.4 画出质控后的小提琴图

plotColData(sce, x="Characteristics.individual.", y="altexps_ERCC_percent", colour_by=I(discard))

plotColData(sce, x="Characteristics.individual.", y="subsets_Mito_percent", colour_by=I(discard))

plotColData(sce, x="Characteristics.individual.", y="sum", colour_by=I(discard))

plotColData(sce, x="Characteristics.individual.", y="detected", colour_by=I(discard))

可见 ERCC 含量过高的都被去除了