smart-seq2 实践

单细胞转录组：smart-seq2

参考资料： Smart-seq2 scRNA小鼠发育学习笔记: //www.greatytc.com/p/23a8fcb12923
服务器：150 :/public/home/kcao/Work/smart-seq2
注意事项：安装最新版本R3.6, 项目中有些包，需要依赖3.6版本。github里面存在多个文件，及其源代码： https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq。详细讲解参考简文：Smart-seq2 scRNA小鼠发育学习笔记: //www.greatytc.com/p/23a8fcb12923
背景：小鼠胚胎细胞6个不同时期的单细胞smart-seq 技术，分析得出发育过程中细胞发育方向，以及标记基因在发育过程中基因基因表达变化。如下图

image.png

1-前言及上游介绍

tips:

10x genomic scRNA-seq 给出R1,R2 两个数据:

R1: 存储的是barcode和UMI信息

R2：真正的测序信息，也就是单端测序
smartseq2 scRNA-seq 给出了R1,R2测序信息：

和普通的RNA-seq类似；可以用hisat2+featureCounts进行操作

1.1 配置conda

# 安装conda
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh 
# 激活
source ~/.bashrc

# 添加镜像
conda config --add channels r
conda config --add channels conda-forge
conda config --add channels bioconda
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

# 创建环境
conda  create -n rna  python=2
conda activate rna
conda install -y  sra-tools trim-galore fastqc multiqc hisat2 samtools subread salmon qualimap

#注销当前的rna环境 
# conda deactivate

下载数据

完整数据sra: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA490198&o=acc_s%3Aa 356.04Gb

下载部分数据测试：

cat SRR_Acc_List.txt
SRR7815790
SRR7815850
SRR7815870
SRR7815890
SRR7815910
SRR7815960
SRR7815980 
SRR7816020 *
SRR7816100
SRR7816120
SRR7816130
SRR7816140 
SRR7816160

加速你的下载

ascp -h 检查效果

# 安装教程
# axel和wget 一样的工具
axel  http://download.asperasoft.com/download/sw/connect/3.7.4/aspera-connect-3.7.4.147727-linux-64.tar.gz

tar zxvf aspera-connect-3.7.4.147727-linux-64.tar.gz

#安装
bash aspera-connect-3.7.4.147727-linux-64.sh
# 然后cd到根目录下看看是不是存在了.aspera文件夹，有的话表示安装成功
cd && ls -a
# 将aspera软件加入环境变量，并激活
echo 'export PATH=~/.aspera/connect/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
# 最后检查ascp是不是能用了
ascp --help

prefetch下载

cat SRR_Acc_List.txt |while read i
do prefetch $i -O `pwd` && echo "** ${i}.sra done **"
done

sra2fq

wkd=/public/home/kcao/Work/smart-seq2/01_raw_data
ls ${wkd}/SRR* | while read id;do \
    (fastq-dump --gzip --split-3 -A `basename $id` -O $wkd/rawdata $id &);done

下载参考基因组注释文件

cd /public/home/kcao/genome_human/mm10/Gencode_data
# 从Gencode下载参考基因组
axel  ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/GRCm38.primary_assembly.genome.fa.gz
# 从Gencode下载参考转录组
axel  ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/gencode.vM23.transcripts.fa.gz
# 下载GTF文件
axel ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/gencode.vM23.annotation.gtf.gz

基于比对的hisat2流程

# 以其中SRR7815790为例
cd /public/home/kcao/Work/smart-seq2/02_mapping
genome_index=/public/home/kcao/genome_human/mm10/Gencode_data/mm10_hista2/genome
fq_dir=/public/home/kcao/Work/smart-seq2/01_raw_data/rawdata
map_dir=/public/home/kcao/Work/smart-seq2/02_mapping

for id  in ${fq_dir}/*_1.fastq.gz;do
    echo " hisat2 -p 10 -x ${genome_index} -1 ${id} -2 ${id/_1.fastq.gz/_2.fastq.gz}  --new-summary -S ${map_dir}/$(basename $id _1.fastq.gz).sam &&
    samtools sort -O bam -@ 10  -o ${map_dir}/$(basename $id _1.fastq.gz).bam ${map_dir}/$(basename $id _1.fastq.gz).sam &&
    samtools index -@ 10 -b ${map_dir}/$(basename $id _1.fastq.gz).bam "|qsub -d ./ -N $(basename $id _1.fastq.gz);done

# 查看比对率
(base) [kcao@login 02_mapping]$ ls *.e[0-9]* |while read id; do echo  -e  "${id}:";grep "Aligned concordantly 1 time" $id| python -c "s=input();print(s.strip())";done
SRR7815790.sra.e2467698:
Aligned concordantly 1 time: 3567726 (79.05%)
SRR7815850.sra.e2467699:
Aligned concordantly 1 time: 4570139 (75.28%)
SRR7815870.sra.e2467700:
Aligned concordantly 1 time: 7540359 (77.94%)
SRR7815890.sra.e2467701:
Aligned concordantly 1 time: 6056489 (80.19%)
SRR7815910.sra.e2467702:
Aligned concordantly 1 time: 4780356 (82.23%)
SRR7815960.sra.e2467703:
Aligned concordantly 1 time: 5449535 (80.69%)
SRR7815980.e2467704:
Aligned concordantly 1 time: 4180821 (74.41%)
SRR7816020.e2467705:
Aligned concordantly 1 time: 18989 (51.17%)
SRR7816100.sra.e2467706:
Aligned concordantly 1 time: 3969194 (75.92%)
SRR7816120.sra.e2467707:
Aligned concordantly 1 time: 4726218 (81.82%)
SRR7816130.sra.e2467708:
Aligned concordantly 1 time: 3698428 (75.50%)
SRR7816140.e2467709:
Aligned concordantly 1 time: 5486285 (79.56%)
SRR7816160.sra.e2467710:
Aligned concordantly 1 time: 6039770 (72.84%)

定量:featurecount

cd /public/home/kcao/genome_human/mm10/Gencode_data
gunzip gencode.vM23.annotation.gtf.gz

gtf=/public/home/kcao/genome_human/mm10/Gencode_data/gencode.vM23.annotation.gtf

cd /public/home/kcao/Work/smart-seq2/03_count
featureCounts -T 10 -p -t exon -g gene_name -a $gtf -o /public/home/kcao/Work/smart-seq2/03_count/hisat2_counts.txt /public/home/kcao/Work/smart-seq2/02_mapping/*.bam 1>featureCounts.log 2>&1

2.根据表达矩阵进行分群

1 首先创造表达矩阵

下载完整的表达矩阵数据：

https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/raw/master/data/female_rpkm.Robj

一会要用的基因列表：https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/data/prot_coding.csv

load(file="female_rpkm.Robj")

## 检测细胞类似数目
> table(sub(pattern="_.*",replacement="",x=colnames(female_rpkm) ))
E10.5 E11.5 E12.5 E13.5 E16.5    P6 
   68   100   104    99    85   108

## 查看重复细胞
#（例如：同一个细胞建库两次，这里作者用“rep”进行了标记） 
grep("rep",colnames(female_rpkm))
colnames(female_rpkm)[256:257]

## 去除重复细胞
female_rpkm <- female_rpkm[,!colnames(female_rpkm) %in% grep("rep",colnames(female_rpkm), value=TRUE)]


## 只保留编码基因(去掉类似：X5430419D17Rik、BC003331等)
prot_coding_genes <- read.csv(file="prot_coding.csv")
head(prot_coding_genes)

females <- female_rpkm[rownames(female_rpkm) %in% as.vector(prot_coding_genes$x),]
save(females,file = 'female_rpkm.Rdata')

2 然后使用包装好的代码进行tSNE

2.1 对细胞操作->细胞发育时期的获取

细胞是从6个时间点取出的，于是先找到这6个时间点

load('/female_rpkm.Rdata')

> dim(females)
[1] 21083   563

> head(colnames(females))
[1] "E10.5_XX_20140505_C01_150331_1" "E10.5_XX_20140505_C02_150331_1"
[3] "E10.5_XX_20140505_C03_150331_1" "E10.5_XX_20140505_C04_150331_2"
[5] "E10.5_XX_20140505_C06_150331_2" "E10.5_XX_20140505_C07_150331_3


## 取下划线分隔的第一部分,查看聚类时候，细胞比例
library(tidyverse)

> sub(pattern="_.*",x=colnames(females) ,replacement="") %>% table()
.
E10.5 E11.5 E12.5 E13.5 E16.5    P6 
   68   100   103    99    85   108

2.2 对基因操作->基因过滤与统计

去掉在所有细胞都不表达的基因

> (dim(females))
[1] 21083   563

> females <- females[!rowSums(females)==0,] 
> dim(females)
[1] 16765   563

可以看到去掉了4000多个

计算各种统计指标

# 利用apply函数对每行（每个基因）进行统计
mean_per_gene <- apply(females, 1, mean, na.rm = TRUE) 
sd_per_gene <- apply(females, 1, sd, na.rm = TRUE) 
mad_per_gene <- apply(females, 1, mad, na.rm = TRUE) 
cv = sd_per_gene/mean_per_gene
library(matrixStats)
var_per_gene <- rowVars(as.matrix(females))
cv2=var_per_gene/mean_per_gene^2
# 存储统计结果
cv_per_gene <- data.frame(mean = mean_per_gene,
                          sd = sd_per_gene,
                          mad=mad_per_gene,
                          var=var_per_gene,
                          cv=cv,
                          cv2=cv2)
rownames(cv_per_gene) <- rownames(females)

> head(cv_per_gene)
              mean           sd          mad          var        cv        cv2
eGFP  1.813189e+03 2753.3822476 1071.2285596 7.581114e+06  1.518530   2.305933
Gnai3 8.390717e+01  127.3289363   42.6915607 1.621266e+04  1.517498   2.302799
Cdc45 1.682541e+01   59.5934190    0.0000000 3.551376e+03  3.541870  12.544840
Scml2 1.003519e+00    6.5980242    0.0000000 4.353392e+01  6.574890  43.229175
Apoh  1.384997e-02    0.1541287    0.0000000 2.375567e-02 11.128451 123.842413
Narf  2.097623e+01   49.8831372    0.1365001 2.488327e+03  2.378079   5.655260

# 根据表达量过滤统计结果
> cv_per_gene=cv_per_gene[cv_per_gene$mean>1,]
> dim(cv_per_gene)
[1] 11398     6


# 简易的可视化
with(cv_per_gene,plot(log10(mean),log10(cv2)))

image.png

sd值，它和均值mean、方差var一样，都是对一维数据进行的分析，如果出现两组数据测量尺度差别太大或数据量纲存在差异的话，直接用标准差就不合适了
CV变异系数就可以解决这个问题，它利用原始数据标准差和原始数据平均值的比值来各自消除尺度与量纲的差异。

复杂一点的统计可视化：

其实就是求每列之间的相关性

library(psych)
pairs.panels(cv_per_gene, 
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = TRUE # show correlation ellipses
)

可以得到不同统计指标的关系

2020年3月19日18:29:58

再用作者包装的函数：`getMostVarGenes()`

https://github.com/IStevant/XX-XY-mouse-gonad-scRNA-seq/blob/master/scripts/analysis_functions.R

getMostVarGenes <- function(
        data=data,              # RPKM matrix
        fitThr=1.5,             # Threshold above the fit to select the HGV
        minMeanForFit=1         # Minimum mean gene expression level
        ){
    # Remove genes expressed in no cells
    data_no0 <- as.matrix(
        data[rowSums(data)>0,]
    )
    # Compute the mean expression of each genes
    meanGeneExp <- rowMeans(data_no0)
    names(meanGeneExp)<- rownames(data_no0)

    # Compute the squared coefficient of variation
    varGenes <- rowVars(data_no0)
    cv2 <- varGenes / meanGeneExp^2
 
    # Select the genes which the mean expression is above the expression threshold minMeanForFit
    useForFit <- meanGeneExp >= minMeanForFit

    # Compute the model of the CV2 as a function of the mean expression using GLMGAM
    fit <- glmgam.fit( cbind( a0 = 1, a1tilde = 1/meanGeneExp[useForFit] ), cv2[useForFit] )
    a0 <- unname( fit$coefficients["a0"] )
    a1 <- unname( fit$coefficients["a1tilde"])
 
    # Get the highly variable gene counts and names
    fit_genes <- names(meanGeneExp[useForFit])
    cv2_fit_genes <- cv2[useForFit]
    fitModel <- fit$fitted.values
    names(fitModel) <- fit_genes
    HVGenes <- fitModel[cv2_fit_genes>fitModel*fitThr]
    print(length(HVGenes))

    # Plot the result
    plot_meanGeneExp <- log10(meanGeneExp)
    plot_cv2 <- log10(cv2)
    plotData <-  data.frame(
        x=plot_meanGeneExp[useForFit],
        y=plot_cv2[useForFit],
        fit=log10(fit$fitted.values),
        HVGenes=log10((fit$fitted.values*fitThr))
    )
    p <- ggplot(plotData, aes(x,y)) +
    geom_point(size=0.1) +
    geom_line(aes(y=fit), color="red") +
    geom_line(aes(y=HVGenes), color="blue") +
    theme_bw() +
    labs(x = "Mean expression (log10)", y="CV2 (log10)")+
    ggtitle(paste(length(HVGenes), " selected genes", sep="")) +
    theme(
        axis.text=element_text(size=16),
        axis.title=element_text(size=16),
        legend.text = element_text(size =16),
        legend.title = element_text(size =16 ,face="bold"),
        legend.position= "none",
        plot.title = element_text(size=18, face="bold", hjust = 0.5),
        aspect.ratio=1
    )+
    scale_color_manual(
        values=c("#595959","#5a9ca9")
    )
    print(p)

    # Return the RPKM matrix containing only the HVG
    HVG <- data_no0[rownames(data_no0) %in% names(HVGenes),]
    return(HVG)
}


females_data <- getMostVarGenes(females, fitThr=2)
> dim(females_data)
[1] 822 563

找到了822个变化比较大的基因

2020年3月19日20:35:55

females_data <- log(females_data+1)
> females_data[1:4,1:4]
         E10.5_XX_20140505_C01_150331_1 E10.5_XX_20140505_C02_150331_1
Ngfr                                  0                              0
Slc22a18                              0                              0
Tspan32                               0                              0
Gmpr                                  0                              0
         E10.5_XX_20140505_C03_150331_1 E10.5_XX_20140505_C04_150331_2
Ngfr                          0.4204863                       3.619946
Slc22a18                      0.0000000                       0.000000
Tspan32                       0.0000000                       0.000000
Gmpr                          0.0000000                       0.000000

save(females_data,file = 'females_hvg_matrix.Rdata')

2.3 6个发育时期RtSNE分析

t-SNE 过程很耗时，先进行PCA 选取变化最显著的主成分，再通过t-SNE 降维可视化。

先是PCA

针对上面的822个HVGs进行操作

female_sub_pca <- FactoMineR::PCA(
  t(females_data), 
  ncp = ncol(females_data), 
  graph=FALSE
)

然后挑选最显著的主成分，作为tSNE的输入,结果选取了9个维度进行T-SNE可视化

significant_pcs <- jackstraw::permutationPA(
  female_sub_pca$ind$coord, 
  B = 100, 
  threshold = 0.05, 
  verbose = TRUE, 
  seed = NULL
)$r
> significant_pcs
[1] 9

然后使用上面`jackstraw`挑出的显著主成分进行tSNE

# 6个时期给定6个颜色
female_stagePalette <-  c(
  "#2754b5", 
  "#8a00b0", 
  "#d20e0f", 
  "#f77f05", 
  "#f9db21",
  "#43f14b"
)
female_t_sne <- run_plot_tSNE(
  pca=female_sub_pca,
  pc=significant_pcs,
  iter=5000,
  conditions=female_stages,
  colours=female_stagePalette
)

2020年3月19日21:00:06

2.4 根据PCA结果进行层次聚类

采用的方法是：Hierarchical Clustering On Principle Components (HCPC)

用PCA降维的9个主成分，降低数据复杂度。再根据经验，检测层次聚类，分成四类。

# 使用9个显著主成分重新跑PCA
res.pca <- FactoMineR::PCA(
  t(females_data), 
  ncp = significant_pcs, 
  graph=FALSE
)
# 作者根据经验认为分成4群比较好解释，于是设置4
res.hcpc <- FactoMineR::HCPC(
  res.pca, 
  graph = FALSE,
  min=4
)

# 得到分群结果
female_clustering <- res.hcpc$data.clust$clust


> table(female_clustering)
female_clustering
  1   2   3   4 
 90 240 190  43 
 
 # 重新命名
female_clustering <- paste("C", female_clustering, sep="")
names(female_clustering) <- rownames(res.hcpc$data.clust)
# 将C1和C2调换位置
female_clustering[female_clustering=="C1"] <- "C11"
female_clustering[female_clustering=="C2"] <- "C22"
female_clustering[female_clustering=="C22"] <- "C1"
female_clustering[female_clustering=="C11"] <- "C2"
> table(female_clustering)
female_clustering
 C1  C2  C3  C4 
240  90 190  43 


write.csv(female_clustering, file="female_clustering.csv")

2020年3月19日21:16:58

3.标记基因可视化

使用包装好的代码，对我们研究的marker基因，通过降维查看marker基因在降维后的空间分布

同时通过marker基因在空间的分布图，可以验证前面的细胞聚类的正确性

3.1 加载表达矩阵、获得cluster信息

####################################
# 3.1 加载表达矩阵、获得cluster信息 
####################################


rm(list = ls()) 
options(warn=-1) 
options(stringsAsFactors = F)
source("../04_cluster/analysis_functions.R")
load('../04_cluster/female_rpkm.Rdata')

# 加载之前HCPC分群结果
cluster <- read.csv('../04_cluster/female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]

>table(female_clustering)
female_clustering
 C1  C2  C3  C4 
240  90 190  43

3.2 拿到文章中的marker基因列表

作者对一些marker基因进行可视化，可以从github上，下载基因集

####################################
# 3.2 拿到文章中的marker基因列表
####################################
# 作者选择的14个marker基因
markerGenes <- c(
  "Nr2f1",
  "Nr2f2",
  "Maf",
  "Foxl2",
  "Rspo1",
  "Lgr5",
  "Bmp2",
  "Runx1",
  "Amhr2",
  "Kitl",
  "Fst",
  "Esr2",
  "Amh",
  "Ptges"
)

3.3 提取marker基因的表达矩阵

矩阵：行是每一个marker基因，列是每一个单细胞名称

####################################
# 3.3 提取marker基因小表达矩阵
####################################
gene_subset <- as.matrix(log(females[rownames(females) %in% markerGenes,]+1))

dim(gene_subset)
gene_subset[1:4,1:4]

cl1_gene_subset <- gene_subset[, colnames(gene_subset) %in% names(female_clustering[female_clustering=="C1"])]
cl2_gene_subset <- gene_subset[, colnames(gene_subset) %in% names(female_clustering[female_clustering=="C2"])]
cl3_gene_subset <- gene_subset[, colnames(gene_subset) %in% names(female_clustering[female_clustering=="C3"])]
cl4_gene_subset <- gene_subset[, colnames(gene_subset) %in% names(female_clustering[female_clustering=="C4"])]


heatmap_gene_subset <- cbind(
  cl1_gene_subset, 
  cl2_gene_subset, 
  cl3_gene_subset,
  cl4_gene_subset
)

3.4 根据marker基因的顺序，重新排列这个矩阵

marker按照文章的顺序进行排列--对行进行排序

####################################
# 3.4 根据marker基因的顺序，重新排列这个矩阵
####################################
# 之前是这样
rownames(heatmap_gene_subset);
markerGenes

# 得到marker基因在heatmap_gene_subset中的位置
match(markerGenes,rownames(heatmap_gene_subset))

# 然后就能提取出和marker基因顺序一样的heatmap_gene_subset
heatmap_gene_subset <- heatmap_gene_subset[match(markerGenes,rownames(heatmap_gene_subset)),]

# 之后是这样
rownames(heatmap_gene_subset);markerGenes

[1] "Nr2f1" "Nr2f2" "Maf"   "Foxl2" "Rspo1" "Lgr5"  "Bmp2"  "Runx1" "Amhr2" "Kitl"  "Fst"   "Esr2"  "Amh"   "Ptges"
 [1] "Nr2f1" "Nr2f2" "Maf"   "Foxl2" "Rspo1" "Lgr5"  "Bmp2"  "Runx1" "Amhr2" "Kitl"  "Fst"   "Esr2"  "Amh"   "Ptges"

3.5 修改表达矩阵的列名，统计每一个时期细胞数目

####################################
### 3.5 修改表达矩阵的列名，得到6个时间点信息
####################################
heatmap_female_stages <- sapply(strsplit(colnames(heatmap_gene_subset), "_"), `[`, 1)
>table(heatmap_female_stages)
heatmap_female_stages
E10.5 E11.5 E12.5 E13.5 E16.5    P6 
   68   100   103    99    85   108

3.6 用包装好的pheatmap画热图

####################################
### 1.6 用包装好的pheatmap画热图
####################################

# 看到H图中，列被分成了4栏，那么这个就是根据colbreaks来划分的。colbreaks的意思就是从哪里到哪里这是一块。当有多个分组又想画一个分割线的话，这个参数就很有用
colbreaks <- c(
  ncol(cl1_gene_subset),
  ncol(cl1_gene_subset)+ncol(cl2_gene_subset), 
  ncol(cl1_gene_subset)+ncol(cl2_gene_subset)+ncol(cl3_gene_subset)
)

# 然后就是上色，6个时间点和4个群使用自定义的颜色
cluster_color <- c(
  C1="#560047",
  C2="#a53bad", 
  C3="#eb6bac", 
  C4="#ffa8a0"
)

stage_color=c(
  E10.5="#2754b5", 
  E11.5="#8a00b0", 
  E12.5="#d20e0f", 
  E13.5="#f77f05", 
  E16.5="#f9db21",
  P6="#43f14b"
)

# 开始画热图
library(pheatmap)
png("female_marker_heatmap.png")
plot_heatmap_2(
  heatmap_gene_subset, 
  female_clustering, 
  heatmap_female_stages, 
  rowbreaks, 
  colbreaks,
  cluster_color,
  stage_color
)
dev.off()

行按照标记基因顺序，列按照细胞类型及其细胞时期排序

image.png

3.7 用包装好的ggboxplot画小提琴图

####################################
### 1.7 用包装好的ggboxplot画小提琴图
####################################
# 小提琴图

pdf("step2.1-B-markers-violin.pdf", width=10, height=22)
require(gridExtra)

# 每个基因的小提琴图都有4个cluster，对它们用不同的颜色
female_clusterPalette <- c(
  "#560047", 
  "#a53bad", 
  "#eb6bac", 
  "#ffa8a0"
)
# 每个基因做一个小提琴图，并用for循环保存在p这个列表中
p <- list()
for (genes in markerGenes) {
  p[[genes]] <- violin_gene_exp(
    genes, 
    females, 
    female_clustering, 
    female_clusterPalette
  )
}
# 最后组合起来，每列显示3张
do.call(grid.arrange,c(p, ncol=3))
dev.off()

在不同时期，部分的标记基因基因表达量小提琴图

image.png

3.8 用包装好的`geom_point`+`geom_density_2d`画等高线图

########################
### 1.8 用包装好的geom_point+geom_density_2d画等高线图
########################
pdf("step2.1-C-markers-tSNE-density.pdf", width=16, height=28)
require(gridExtra)

load('../04_cluster/female_tsne.Rdata')
female_tsne <-female_t_sne_new_clusters
p <- list()
for (genes in markerGenes) {
  p[[genes]] <- tsne_gene_exp(
    female_tsne,
    genes, 
    females
  )
}

do.call(grid.arrange,c(p, ncol=4))
dev.off()

坐标点为细胞T-sne 降维后位点，颜色深浅代表着表达量高低，等高线代表着细胞分布密度

image.png

4.差异基因分析

前面得到的6个发育时期和4个分群，而且还可视化了一些marker基因，那么现在就要对这4群细胞进行差异分析

使用monocle V2进行差异分析

4.1 准备表达矩阵和分群信息

## 加载包装好的R代码
rm(list = ls()) 
options(warn=-1) 
options(stringsAsFactors = F)
source("../04_cluster/analysis_functions.R")

# 差异分析一般都要求使用count矩阵
load('../04_cluster/female_count.Rdata')

# 6个发育时期获取
head(colnames(female_count))
female_stages <- sapply(strsplit(colnames(female_count), "_"), `[`, 1)
names(female_stages) <- colnames(female_count)

table(female_stages)
# female_stages
# E10.5 E11.5 E12.5 E13.5 E16.5    P6 
# 68   100   103    99    85   108

# 4个cluster获取
cluster <- read.csv('../04_cluster/female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
table(female_clustering)

4.2 利用monocle V2构建对象

需要三样东西：表达矩阵、细胞信息、基因信息，类似进行DESeq2 分析步骤


############################
## 1.2 准备表达矩阵和分群信息
############################

## 表达矩阵
dim(female_count)
## 细胞分群信息（包括6个stage和4个cluster）
table(female_stages)
table(female_clustering)
## 基因信息
gene_annotation <- as.data.frame(rownames(female_count))


################
## 构建分组信息
################
library(monocle) 
# 直接使用作者包装的函数，代码更简洁
DE_female <- prepare_for_DE (
  female_count, 
  female_clustering, 
  female_stages
)

4.3 进行差异分析

使用包装好的函数：differentialGeneTest函数

得到4435个差异基因，并且将这些差异基因分配到4个cluster 中。

################################
# 1.3 进行差异分析
################################
start_time <- Sys.time()
female_DE_genes <- findDEgenes(
  DE_female, 
  qvalue=0.05
)
end_time <- Sys.time()
end_time - start_time

# [1] "4435 significantly DE genes (FDR<0.05)."
# [1] "3268 significantly DE genes (FDR<0.01)."


############
# 定位差异基因所属于细胞cluster
###########
# 有了差异基因以后，继续使用RPKM值
load('../04_cluster/female_rpkm.Rdata')
de_clusters <- get_up_reg_clusters(
  females, 
  female_clustering, 
  female_DE_genes
)

View(de_clusters)

save(de_clusters,file = "step3.1-DEG-monocle_summary.Rdata")

最后一列显示这个差异基因，是哪个cluster 的marker基因

image.png

5.差异基因功能注释

找差异基因的方法，有很多，例如DESeq2,EdgeR,Wilcox,,SCDE 等等，

得到差异基因，及其所属类别，就可以对每一类进行富集分析

5.1 对monocle的结果进行注释,分组转换ID


rm(list=ls())
options(stringsAsFactors = F)
library(clusterProfiler)

load(file = 'step3.1-DEG-monocle_summary.Rdata')
# 得到差异基因名
de_genes <- subset(de_clusters, qval<0.05)
length(de_genes$genes)

# 然后基因ID转换
entrez_genes <- bitr(de_genes$genes, fromType="SYMBOL", 
                     toType="ENTREZID", 
                     OrgDb="org.Mm.eg.db")
dim(entrez_genes)
# 作者剔除掉一个基因名
entrez_genes <- entrez_genes[!entrez_genes$ENTREZID %in% "101055843",]
length(entrez_genes$SYMBOL)

de_gene_clusters :记录了差异基因所属类别。

################################
# 把存在ENTREZID的那些基因的基因名和cluster信息提取出来
de_gene_clusters <- de_genes[de_genes$genes %in% entrez_genes$SYMBOL,
                             c("genes", "cluster")]



> head(de_gene_clusters)
        genes cluster
A4galt A4galt      C2
Aaas     Aaas      C1
Aacs     Aacs      C3
Aagab   Aagab      C3
Aar2     Aar2      C3
Aard     Aard      C4

################################
# 保持de_gene_clusters$genes的顺序不变，将他的symbol变成entrez ID
de_gene_clusters <- data.frame(
  ENTREZID=entrez_genes$ENTREZID[entrez_genes$SYMBOL %in% de_gene_clusters$genes],
  cluster=de_gene_clusters$cluster
)

5.2 对每个cluster的差异基因，进行GO富集分析

用clusterProfiler 对每一个cluster的差异基因，进行可视化。

#######
# 将差异基因对应到每个cluster，拆分成list
list_de_gene_clusters <- split(de_gene_clusters$ENTREZID, 
                               de_gene_clusters$cluster)

#######
# 整理数据结构list， 进行GO 富集
formula_res <- compareCluster(list_de_gene_clusters,
                     fun="enrichGO",
                     OrgDb="org.Mm.eg.db", 
                     ont          = "BP",
                     pAdjustMethod = "BH",
                     pvalueCutoff=0.01,
                     qvalueCutoff  = 0.05)

# 可视化
pdf('DEG_GO_each_cluster.pdf',width = 11,height = 6)
dotplot(formula_res, showCategory=5)
dev.off()

image.png

6.发育谱系推断及可视化

谱系研究适用于研究在细胞分化过程中，细胞发育趋势，同时可以看到标记基因，在这个发育过程中，基因表达会发生显著变化。

所需的数据材料：

表达矩阵，高变化的基因cv，发育时期的信息，分群信息，就可以进行谱系推断。有一些方法可以构建发育谱系，比如DiffusionMap、Slingshot

构建发育谱系，先走一下DiffusionMap的流程，得到几个重要的DC；接着走一下Slingshot函数，就会得到谱系结果

6.1 运行DiffusionMap，进行降维

准备表达矩阵、HVGs、分群信息

source("../04_cluster/analysis_functions.R")

### 1.1 准备表达矩阵、HVGs、分群信息
# 表达矩阵
load('../04_cluster/female_rpkm.Rdata')
dim(females)

# HVGs
load('../04_cluster/females_hvg_matrix.Rdata')
dim(females_data)

# 6个发育时期获取
head(colnames(females))
female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
names(female_stages) <- colnames(females)
table(female_stages)

# 4个cluster获取
cluster <- read.csv('../04_cluster//female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
table(female_clustering)

运行DiffusionMap: 类似PCA一样，找到前几位DC

### 运行
female_dm <- run_diffMap(
  females_data, 
  female_clustering,
  sigma=15
)

### 作图探索
plot_eigenVal(
  dm=female_dm
)

可以看到前3个DC成分，贡献很大

image.png

三个DC所组成的空间，主要就是选取了前3个DC成分，做了三维空间的映射，然后把点的颜色分别按照cluster和stage两种不同的属性上色，
- 下图为4个cluster 空间位置
image.png
- 下图为6个stage 空间位置

image.png

6.2 进行Slingshot分析，计算具体的分化路径

利用DiffusionMap，DC分析结果，认为DC4 为拐点。

image.png

计算最优的分化路线，及其分化方向


dm=female_dm
dim=c(1:4)
condition=factor(female_clustering)

data <- data.frame(
  dm@eigenvectors[,dim]
)

female_lineage <- slingshot(
  data, 
  condition, 
  start.clus = "C1", 
  end.clus=c("C2", "C4"),
  maxit=100000,
  shrink.method="cosine"
  # shrink.method="tricube"
)
# 看下结果
> female_lineage
class: SlingshotDataSet 

 Samples Dimensions
     563          4

lineages: 2 
Lineage1: C1  C3  C4  
Lineage2: C1  C2  

curves: 2 
Curve1: Length: 1.3739  Samples: 453.59
Curve2: Length: 0.74646 Samples: 312.73

可以看出4维情况，推荐的分化线路2条，及其给出了对应的细胞数目

输出落在分化路线上的细胞

# 推断的细胞发育谱系结果
female_pseudotime <- get_pseudotime(female_lineage, wthres=0.9)
rownames(female_pseudotime) <- colnames(females)

结果可以看到数据框中，存在两列，curve1 和curve2.

当curve2 为NA ，表示此细胞只存在于curve1中，也就是分支上，

当curve2 && curve1 !=NA ,表明此细胞落在主枝上。

image.png

进一步用图形解释：图中显示了细胞发育状态及其方向
- 最初是由E10.5作为起点发育的，然后分化成两个方向
- 看到P6这个时期是在两条轨迹中都存在的，说明这个时期的细胞存在异质性。也就是说，虽然都是P6时期取的细胞，但是它们实际的发育方向是不用的
- 这个算法就是帮助我们认识到不同时期内部的各个细胞，它们依然还存在着差异

image.png

6.3 谱系发育相关基因可视化

最初我们知道细胞有6个时期（就是取样的6个时间点）；然后进行聚类发现这些细胞能分成4个cluster（意思就是虽然是一个时间点取的细胞，依然可能属于不同类型）；后来进行谱系推断，又增加了一个细胞属性（就是2条不同的发育轨迹）

载入之前的结果

options(warn=-1) 
options(stringsAsFactors = F)


load('../04_cluster/female_rpkm.Rdata')

对谱系推断结果进行归一化

要让两条谱系可以进行比较，需要对每条轨迹进行标准化，除以个值最大值。

## 第一条
pseudotime_lin <- female_pseudotime[,"curve1"]
max_pseudotime <- max(pseudotime_lin, na.rm = TRUE)
pseudotime_lin1_percent <- (pseudotime_lin*100)/max_pseudotime

## 第二条
pseudotime_lin <- female_pseudotime[,"curve2"]
max_pseudotime <- max(pseudotime_lin, na.rm = TRUE)
pseudotime_lin2_percent <- (pseudotime_lin*100)/max_pseudotime

# 现在female_pseudotime中的两条轨迹结果都在0-100之间了
female_pseudotime[,"curve1"] <- pseudotime_lin1_percent
female_pseudotime[,"curve2"] <- pseudotime_lin2_percent

可视化

由于发育过程是连续的，特异基因会发生明显的表达变化

1.不区分分化路径，只考虑时间

## 先给一个颜色, 6个stage颜色
female_stagePalette <- c(
  E10.5="#2754b5", 
  E11.5="#8a00b0", 
  E12.5="#d20e0f", 
  E13.5="#f77f05", 
  E16.5="#f9db21",
  P6="#43f14b"
)
# 4个cluster颜色
female_clusterPalette <- c(
  C1="#560047",
  C2="#a53bad", 
  C3="#eb6bac", 
  C4="#ffa8a0"
)
# 2个发育谱系颜色
female_clusterPalette2 <- c(
  "#ff6663", 
  "#3b3561"
)
pdf("Amhr2_develop.pdf")
## 做第一个谱系的一个基因（以Amhr2基因为例）
plot_smoothed_gene_per_lineage(
  rpkm_matrix=females, # RPKM表达矩阵
  pseudotime=female_pseudotime,  #谱系推断结果
  lin=c(1), # 对第一个谱系操作
  gene="Amhr2",  #画Amhr2基因变化
  stages=female_stages, # 发育时间点分类
  clusters=female_clustering, # cluster分类
  stage_colors=female_stagePalette,
  cluster_colors=female_clusterPalette,
  lineage_colors=female_clusterPalette2
)
dev.off()

Amhr2 【markerGenes】发育过程表达变化

image.png

区分发育方向，查看两条发育方向，标记基因表达变化情况。以Amhr2 为例：

# 做某个基因在两个谱系中的变化（以Amhr2基因为例）
pdf("Amhr2_develop_exp_two_lineage.pdf")
plot_smoothed_gene_per_lineage(
  rpkm_matrix=females, 
  pseudotime=female_pseudotime, 
  lin=c(1,2),
  gene="Amhr2", 
  stages=female_stages, 
  clusters=female_clustering, 
  stage_colors=female_stagePalette,
  cluster_colors=female_clusterPalette,
  lineage_colors=female_clusterPalette2
)
dev.off()

标记基因Amhr2【markerGene】在两条分化方向上，基因表达变化情况。

看到这个Amhr2基因在第一个谱系中变化很大，尤其是到后期；而在第二个谱系中基本保持平衡，这就说明这个基因就是第一个谱系中重要的基因

image.png

对多个标记基因，考虑分支/不考虑分支进行画图

gene_list <- c(
  "Sall4",
  "Sox11",
  "Gata4",
  "Lgr5",
  "Runx1",
  "Foxl2",
  "Hey2",
  "Wnt5a",
  "Pdgfra",
  "Nr2f2"
)

plot_smoothed_genes <- function(genes, lin){
  female_clusterPalette2 <- c("#ff6663", "#3b3561")
  for (gene in genes){
    plot_smoothed_gene_per_lineage(
      rpkm_matrix=females, 
      pseudotime=female_pseudotime, 
      lin=lin,
      gene=gene, 
      stages=female_stages, 
      clusters=female_clustering, 
      stage_colors=female_stagePalette,
      cluster_colors=female_clusterPalette,
      lineage_colors=female_clusterPalette2
    )
  }
}

pdf("interesting_genes_in_lineage.pdf", width=4, height=4)
plot_smoothed_genes(gene_list, 1) # plot only lineage 1
plot_smoothed_genes(gene_list, 2) # plot only lineage 2
plot_smoothed_genes(gene_list, c(1,2)) # plot the two moleages in the same graph to see the divergence
dev.off()

最后编辑于：2020.03.21 09:05:19

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 212,454评论 6赞 493
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 90,553评论 3赞 385
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 157,921评论 0赞 348
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 56,648评论 1赞 284
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 65,770评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 49,950评论 1赞 291
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,090评论 3赞 410
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,817评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,275评论 1赞 303
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,592评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,724评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,409评论 4赞 333
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,052评论 3赞 316
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,815评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,043评论 1赞 266
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,503评论 2赞 361
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,627评论 2赞 350

smart-seq2 实践

单细胞转录组：smart-seq2

1-前言及上游介绍

1.1 配置conda

下载数据

prefetch下载

sra2fq

下载参考基因组注释文件

基于比对的hisat2流程

定量:featurecount

2.根据表达矩阵进行分群

1 首先创造表达矩阵

2 然后使用包装好的代码进行tSNE

2.1 对细胞操作->细胞发育时期的获取

2.2 对基因操作->基因过滤与统计

去掉在所有细胞都不表达的基因

计算各种统计指标

复杂一点的统计可视化：

再用作者包装的函数：getMostVarGenes()

2.3 6个发育时期RtSNE分析

先是PCA

然后使用上面jackstraw挑出的显著主成分进行tSNE

2.4 根据PCA结果进行层次聚类

3.标记基因可视化

3.1 加载表达矩阵、获得cluster信息

3.2 拿到文章中的marker基因列表

3.3 提取marker基因的表达矩阵

3.4 根据marker基因的顺序，重新排列这个矩阵

3.5 修改表达矩阵的列名，统计每一个时期细胞数目

3.6 用包装好的pheatmap画热图

3.7 用包装好的ggboxplot画小提琴图

3.8 用包装好的geom_point+geom_density_2d画等高线图

4.差异基因分析

使用monocle V2进行差异分析

4.1 准备表达矩阵和分群信息

4.2 利用monocle V2构建对象

需要三样东西：表达矩阵、细胞信息、基因信息，类似进行DESeq2 分析步骤

4.3 进行差异分析

5.差异基因功能注释

5.1 对monocle的结果进行注释,分组转换ID

5.2 对每个cluster的差异基因，进行GO富集分析

6.发育谱系推断及可视化

6.1 运行DiffusionMap，进行降维

准备表达矩阵、HVGs、分群信息

运行DiffusionMap: 类似PCA一样，找到前几位DC

6.2 进行Slingshot分析，计算具体的分化路径

6.3 谱系发育相关基因可视化

载入之前的结果

对谱系推断结果进行归一化

可视化

区分发育方向，查看两条发育方向，标记基因表达变化情况。以Amhr2 为例：

对多个标记基因，考虑分支/不考虑分支进行画图

推荐阅读更多精彩内容

再用作者包装的函数：`getMostVarGenes()`

然后使用上面`jackstraw`挑出的显著主成分进行tSNE

3.8 用包装好的`geom_point`+`geom_density_2d`画等高线图