一、建立参考基因组索引
1、下载参考基因组和注释文件
2、构建适用于cellranger分析所用的参考基因组和注释文件索引（index），代码如下：
##根据自己分析的物种（此处分析的是猪的数据，其他物种同理）的注释文件提取gtf信息
cellranger mkgtf chr_susScr11.gtf chr_susScr11.filtered.gtf --attribute=gene_biotype:protein_coding
##根据自己分析的物种参考基因组和注释文件构建gtf索引
cellranger mkref --genome=chr_susScr11_cellranger \
--nthreads=10 \
--fasta=chr_susScr11.fa \
--genes=chr_susScr11.filtered.gtf
二、比对和计数分析
1、输入文件
Cellranger count 的fastq输入文件格式必须要求https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/fastq-input
格式如下：
[Sample Name]S1_L00[Lane Number][Read Type]_001.fastq.gz
Where Read Type is one of:
I1: Sample index read (optional)*
R1: Read 1
R2: Read 2
注意：1）cellranger count的输入文件格式必须满足上述要求；2）后缀必须是fastq或者fastq.gz；如下
image-20200502201843961.png
2、构建shell流程批量分析
2 ##Pig single cell data analysis
3 ##########QC 1###########
4 genomedir=/data/youhua_data/Database/susScr11 ##参考基因组index目录
5 datadir=/data/youhua_data/Data/WRF_singlecell ##文件目录
6
7
8 sample='Fat1 Fat2 Fat3 Fat4' ##sample名
9 date
10 for s in $sample
11 do
12 date
13 cellranger count --id=${s}_cellrange_out \ ##cellranger count计数分析
14 --fastqs=$datadir/trim_output_dir/ \
15 --sample=$s \
16 --transcriptome=$genomedir/chr_susScr11_cellranger
17 date
18 wait
19 done
20 exit
3、输出文件
1)输出文件夹目录
image-20200502202226618.png
2）outs目录文件
image-20200502202350207.png
三、Seurat分析
library(dplyr)
library(Seurat)
pig.data <- Read10X(data.dir = "outs/filtered_feature_bc_matrix/") ##加载cellranger cout分析后的数据
##创建对象，min.feature为基因；counts为rawdata（10x genomic）或者TPM
pig <- CreateSeuratObject(counts = pig.data, project = "pig_singlecell", min.cells = 3, min.features = 200)
head(pig)
pig
pig[["percent.mt"]] <- PercentageFeatureSet(pig, pattern = "^MT-") ##添加线粒体基因表达属性
png("vlnplot.png",width=1000,height=600)
##绘制单个细胞的基因表达个数、基因表达量以及线粒体基因百分比之间的关系，用于筛选过滤细胞（死细胞等离群细胞）
VlnPlot(pig, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) ##nFeature_RNA为单个细胞的基因表达个数，nCount_RNA基因表达量，percent.mt为线粒体基因百分比；
dev.off()
##绘制nCount_RNA、percent.mt以及nFeature_RNA的相关性散点图，用于筛选过滤细胞
png("boxplot.png",width=1000,height=600)
plot1 <- FeatureScatter(pig, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pig, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
dev.off()
##根据上面分析结果确定筛选条件
pig <- subset(pig, subset = nFeature_RNA > 200 & nFeature_RNA < 1500 & percent.mt < 5) ##过滤：基因数目大于200，不大于2500，MT-percent
##取Log进行数据标准化处理，
pig <- NormalizeData(pig, normalization.method = "LogNormalize", scale.factor = 10000)
##寻找变异系数较大的基因进行分析
pig <- FindVariableFeatures(pig, selection.method = "vst", nfeatures = 1500)
##可视化变异系数较大的细胞
png("feature.png",width=1000,height=600)
par(mfrow=c(1,2))
top10 <- head(VariableFeatures(pig), 10)
plot1 <- VariableFeaturePlot(pig)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = FALSE)
CombinePlots(plots = list(plot1, plot2))
dev.off()
##
all.genes <- rownames(pig)
pig <- ScaleData(pig, features = all.genes)
pig <- RunPCA(pig, features = VariableFeatures(object = pig))
print(pig[["pca"]], dims = 1:5, nfeatures = 5)
png("dim.png")
VizDimLoadings(pig, dims = 1:2, reduction = "pca")
DimPlot(pig, reduction = "pca")
dev.off()
png("dimheatmap.png")
DimHeatmap(pig, dims = 1:9, cells = 500, balanced = TRUE)
dev.off()
pig <- JackStraw(pig, num.replicate = 100)
pig <- ScoreJackStraw(pig, dims = 1:20)
png("pc_jack.png")
JackStrawPlot(pig, dims = 1:20)
dev.off()
png("elbow.png")
ElbowPlot(pig)
dev.off()
pig <- FindNeighbors(pig, dims = 1:15)
pig <- FindClusters(pig, resolution = 0.4)
pig <- RunUMAP(pig, dims = 1:15)
png("umap.png")
DimPlot(pig, reduction = "umap")
dev.off()
saveRDS(pig, file = "./pig_tutorial.rds")
cluster1.markers <- FindMarkers(pig, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
cluster5.markers <- FindMarkers(pig, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
pig.markers <- FindAllMarkers(pig, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pig.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
write.table(file='tsne.txt' ,Embeddings(pig@reductions$umap),sep="\t",col.names=F,quote=F)
write.table(file='tsne.class' ,pig@active.ident,sep="\t",col.names=F,quote=F)
write.table(pig.markers,sep="\t",file="markers.xlsx",row.names=F)
cluster1.markers <- FindMarkers(pig, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
png("de1.png")
VlnPlot(pig, features = c("DES", "HSPB7"))
dev.off()
png("de2.png")
VlnPlot(pig, features = c("UBE2C", "CKS2"), slot = "counts", log = TRUE)
dev.off()
#FeaturePlot(pig, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
png("gene.png")
FeaturePlot(pig, features = c("ACTA2", "TAGLN", "THBS1", "TPM2", "DES", "CRYAB", "TAGLN", "MYL9","CNN1"))
dev.off()
png("heat.png")
top10 <- pig.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pig, features = top10$gene) + NoLegend()
dev.off()
