关于如何用二代测序数据注释单细胞细胞群

cleandata=/home/project/data/cleandata/trim_galore
rawdata=/home/project/data/rawdata/
cat sample.ID | while read id
do
  trim_galore --phred33 -q 20  --stringency 3 --fastqc --paired --max_n 3 -o ${cleandata} ${rawdata}/${id}_1.fastq.gz ${rawdata}/${id}_2.fastq.gz
done

index=refdir/server/reference/index/hisat2.2.1/hg38/genome
inputdir=/home/project/data/cleandata/trim_galore
outdir=/home/project/data/Mapping/Hisat2
cat /teach/t_rna/data/airway/sra/sample.ID | while read id
do
  hisat2 -p 5 -x ${index} -1 ${inputdir}/${id}_1_val_1.fq.gz -2 ${inputdir}/${id}_2_val_2.fq.gz 2>${id}.log  | samtools sort -@ 3 -o ${outdir}/${id}.Hisat_aln.sorted.bam -  &&  samtools index ${outdir}/${id}.Hisat_aln.sorted.bam ${outdir}/${id}.Hisat_aln.sorted.bam.bai
done

gtf=/home/data/refdir/server/reference/gtf/ensembl/Homo_sapiens.GRCh38.98.gtf.gz
inputdir=/home/project/data/Mapping/Hisat2
# featureCounts对bam文件进行计数
featureCounts -T 6 -p -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.sorted.bam
# 得到表达矩阵
cat all.id.txt | cut -f1,7- > counts.txt
# 处理表头部分省略，大家自己根据需要使用sed处理即可

library(stringr)
library(Seurat)
#1.首先处理单细胞数据，获得各细胞群的平均表达矩阵
{
    #以下所用的文件去文章对应的GSE144320下载即可，作者的表达矩阵颠倒了
    counts=read.csv(file = "input/GSM4285803_scRNA_RawCounts.csv",row.names = 1)
    counts=as.matrix(counts)
    counts=t(counts)
    counts=as.data.frame(counts)
    meta=read.csv(file="input/GSM4285803_scRNA_metaInfo.csv",row.names = 1)
    identical(colnames(counts),rownames(meta))
    #获得了表达矩阵和患者信息，构建seurat对象
    sce <- CreateSeuratObject(counts, meta.data = meta)
    sce #21888 features across 4614 samples within 1 assay 
    #查看sce后，发现这是作者已经过滤后的数据
    sce <- NormalizeData(sce)
    sce <- FindVariableFeatures(sce, nfeatures = 3000)
    sce = ScaleData(sce) 
    sce <- RunPCA(sce, npcs = 30)
    sce <- RunTSNE(sce, dims = 1:30)
    sce <- RunUMAP(sce, dims = 1:30)
    sce <- FindNeighbors(sce, dims = 1:30)
    #这里resolution取0.2才能获得和作者一样的细胞群数目
    sce <- FindClusters(sce,  resolution = 0.2, algorithm = 1)
    DimPlot(sce)
    #最终获得了平均值的表达矩阵，将作者定义好的细胞群，和我们自己分好的细胞群进行保存
    av <-AverageExpression(sce, assays = "RNA",group.by = "CellType") 
    av2 = AverageExpression(sce, assays = "RNA",group.by = "RNA_snn_res.0.2") 
    av=av[["RNA"]];av2=av2[["RNA"]]
}

#2.处理作者使用的CCL和GEO上的公共数据表达矩阵
{
    rawcount1=read.table("input/rawcounts.txt",row.names = 1,sep = "\t",header = T)
    rawcount2=read.table("input/GSM1867011_HEK_293_C1_USC_RNA_40.csv",row.names = 1,sep = "\t",header = T)
    rawcount3=read.table("input/GSM1867012_HEK_293_C2_USC_RNA_41.csv",row.names = 1,sep = "\t",header = T)
    identical(rownames(rawcount2),rownames(rawcount3))
    rawcount23=cbind(rawcount2,rawcount3)
    colnames(rawcount23)=c("GSM1867011","GSM1867012")
    rownames(rawcount23)=str_sub(rownames(rawcount23),1,15)#这个需要提前确定，.号前面的ENSG编号没有重复
    rawcount23=as.data.frame(rowSums(rawcount23)/2)
    colnames(rawcount23)="HEK293T";colnames(rawcount1)=c("MCF10A","SUM44","MDAMB134","MCF7","T47D")
    #以上是rawcount矩阵，现在转化为cpm表达矩阵
    #扫了一遍发现，HEK293T的rawcount下载下来后数值都很高，而且没有0值，猜测作者在上传前就已经对原始矩阵进行过过滤
    #既然如此，我们就不对我们的表达矩阵进行过滤，与HEK293T的rawcount进行合并即可
    library(edgeR)
    tmp=intersect(rownames(rawcount1),rownames(rawcount23))
    length(tmp);dim(rawcount1);dim(rawcount23)
    #获得最后的合并后的表达矩阵，进行cpm转换，需要log，不然差异过大，见https://cloud.tencent.com/developer/article/1625229
    rawcount=cbind(rawcount1[tmp,],rawcount23[tmp,])
    colnames(rawcount)[6]="HEK293T"
    cpmcount=log2(cpm(rawcount)+1)
    #进行基因id的转换
    library(clusterProfiler)
    library(org.Hs.eg.db)
    id2symbol <- bitr(rownames(cpmcount), 
                      fromType = "ENSEMBL", 
                      toType = "SYMBOL", 
                      OrgDb = org.Hs.eg.db)
    dim(id2symbol)
    table(duplicated(id2symbol$SYMBOL))#发现有重复,去重
    id2symbol = id2symbol[!duplicated(id2symbol$SYMBOL),]
    dim(id2symbol)
    cpmcount=cpmcount[id2symbol$ENSEMBL,]
    rownames(cpmcount)=id2symbol$SYMBOL
    save(cpmcount,av,sce,av2,file="Rdata/section8.Rdata")
}

{
    #首先是作者已经定义好的细胞群
    load(file="Rdata/section8.Rdata")
    head(av)
    av=as.data.frame(av);cpmcount=as.data.frame(cpmcount)
    dim(av);dim(cpmcount)
    av=av[intersect(rownames(av),rownames(cpmcount)),]
    cpmcount=cpmcount[intersect(rownames(av),rownames(cpmcount)),]
    dim(av);dim(cpmcount)
    data_merged=cbind(av,cpmcount)
    tt <- cor(data_merged,method = "pearson")
    tt1 <- tt[1:8,9:14]
    pheatmap::pheatmap(tt1) 
    #接着是我们自己未定义过的细胞群
    head(av2)
    av2=as.data.frame(av2);cpmcount=as.data.frame(cpmcount)
    dim(av2);dim(cpmcount)
    av2=av2[intersect(rownames(av2),rownames(cpmcount)),]
    cpmcount=cpmcount[intersect(rownames(av2),rownames(cpmcount)),]
    dim(av2);dim(cpmcount)
    data_merged=cbind(av2,cpmcount)
    #进行相关性分析
    tt <- cor(data_merged,method = "pearson")
    tt1 <- tt[1:8,9:14]
    pheatmap::pheatmap(tt1,cluster_rows = F,cluster_cols = F)
}

pheatmap::pheatmap(tt1,cluster_rows = F,cluster_cols = F,scale = "column") 
#这里我确实不知道该怎么解释，请大家自行理解下