Smartseq2 scRNA小鼠發(fā)育學(xué)習(xí)筆記-1-前言及上游介紹

劉小澤寫于19.10.23
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索Smartseq2技術(shù)及發(fā)育相關(guān)的分析
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
這次會介紹文章的兩張主圖咆槽，還有上游分析尽纽，算是引言部分炊汤。對應(yīng)視頻第三單元1-4講

前言

這次要重復(fù)的文章是：Dissecting Cell Lineage Specification and Sex Fate Determination in Gonadal Somatic Cells Using Single-Cell Transcriptomics

文章方法思路

  • 為了研究卵巢發(fā)育過程中體細(xì)胞譜系特征三妈，作者結(jié)合Fluidigm C1分選和smartseq2建庫史煎，在Tg(Nr5a1-GFP)轉(zhuǎn)基因小鼠性腺分化的6個(gè)時(shí)期（E10.5, E11.5, E12.5, E13.5, E16.5, and post-natalday6）提取了663個(gè)細(xì)胞（GSE119766）阔籽，然后使用Hiseq2000 進(jìn)行雙端100bp測序吱肌，每個(gè)細(xì)胞測10M reads。

  • 經(jīng)過比對卷拘、定量得到表達(dá)矩陣喊废，然后對細(xì)胞進(jìn)行過濾，最后剩下563個(gè)栗弟。進(jìn)行PCA污筷，用Jackstraw提取了最顯著的幾個(gè)主成分，并用FactoMineR進(jìn)行HCPC 聚類乍赫，使用Rtsne 進(jìn)行tSNE瓣蛀，然后分群，得到了4個(gè)胚胎發(fā)育時(shí)期（C1到C4）耿焊，之后對標(biāo)記基因進(jìn)行可視化（包括等高線圖和熱圖揪惦，使用了自己包裝的大量代碼）。

  • 然后對C1-C4每個(gè)群進(jìn)行差異分析（使用monocle2）罗侯，差異基因進(jìn)行功能注釋（使用clusterProfiler和PANTHER）

  • 根據(jù)高變化基因進(jìn)行細(xì)胞譜系推斷，利用diffusion map的Destiny 和Slingshot 溪猿。后來譜系發(fā)育基因進(jìn)行了可視化钩杰、聚類和注釋

綜上，文章得到了6個(gè)發(fā)育時(shí)期诊县、4群細(xì)胞讲弄、2個(gè)發(fā)育軌跡這三種細(xì)胞屬性

重要的信息如下：

分群信息（對應(yīng)課程下游分析之第5、6講）
分群信息
標(biāo)記基因可視化（對應(yīng)課程下游分析之第7依痊、8講）
標(biāo)記基因可視化
差異分析及注釋（對應(yīng)課程下游分析之第9避除、10講）
差異分析及注釋
比較不同的差異分析方法，如monocle胸嘁、ROST（對應(yīng)課程下游分析之第11講）
譜系推斷（對應(yīng)課程下游分析之第12講）
譜系推斷
譜系發(fā)育基因可視化（對應(yīng)課程下游分析之第13講）
譜系發(fā)育基因可視化
譜系發(fā)育基因聚類和注釋（對應(yīng)課程下游分析之第14講）
譜系發(fā)育基因聚類和注釋

上游分析

將會對應(yīng)課程的全部上游分析第1-4講
這部分和轉(zhuǎn)錄組分析很相似瓶摆，所以也是簡單帶過

smartseq2得到的文件和10X的不一樣，10X的數(shù)據(jù)雖然也有R1性宏、R2兩個(gè)數(shù)據(jù)群井，但第一個(gè)存儲的是barcode和UMI信息，而只有第二個(gè)才是真正的測序信息毫胜，也就是單端測序书斜；而smartseq2得到的兩個(gè)R1诬辈、R2都是測序信息，于是它的操作和我們常規(guī)bulk轉(zhuǎn)錄組是類似的荐吉”涸悖可以用hisat2+featureCounts進(jìn)行操作；或者如果不想比對的話样屠，也可以用salmon進(jìn)行操作

配置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

# 創(chuàng)建環(huán)境
conda  create -n rna  python=2
conda activate rna
conda install -y  sra-tools trim-galore fastqc multiqc hisat2 samtools subread salmon qualimap

#注銷當(dāng)前的rna環(huán)境 
# conda deactivate

下載數(shù)據(jù)

XX全部的數(shù)據(jù)比較大穿撮，數(shù)據(jù)在：https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA490198&o=acc_s%3Aa

因此簡單下幾個(gè)體驗(yàn)一下就好，我這里選取了這些：

# wkd=/YOUR_PATH/
# cd $wkd/sra
# 這些數(shù)據(jù)全部都放在了EBI瞧哟，而不是NCBI
SRR7815790
SRR7815850
SRR7815870
SRR7815890
SRR7815910
SRR7815960
SRR7815980
SRR7816020
SRR7816100
SRR7816120
SRR7816130
SRR7816140
SRR7816160

如果prefetch + ascp可以用的話混巧，那是最好，因?yàn)槟鞘亲钍∈碌南螺d辦法勤揩，而且不用管SRA數(shù)據(jù)是來自歐洲的EBI還是美國的NCBI咧党。如果出現(xiàn)了使用https下載方式，而非快速的fasp下載陨亡，可以看這個(gè)：來吧傍衡，加速你的下載

需要注意EBI和NCBI的SRA數(shù)據(jù)存儲結(jié)構(gòu)有些區(qū)別：

  • NCBI的文件存在ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/

  • EBI文件存在ftp://ftp.sra.ebi.ac.uk/vol1/srr

    例如要下載SRR7815792.sra，它就是ftp://ftp.sra.ebi.ac.uk/vol1/srr/SRR781/002/SRR7815792负蠕，倒數(shù)第二個(gè)位置的002是和SRR ID的最后一個(gè)數(shù)字對應(yīng)的蛙埂，另外EBI下載的名稱沒有后綴.sra

sra2fq

wkd=/YOUR_PATH/
ls $wkd/sra/SRR* | while read i;do \
    (fastq-dump --gzip --split-3 -A `basename $i` -O $wkd/rawdata $i &);done

下載參考數(shù)據(jù)

# 從Gencode下載參考基因組
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/GRCm38.primary_assembly.genome.fa.gz
# 從Gencode下載參考轉(zhuǎn)錄組
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/gencode.vM23.transcripts.fa.gz
# 下載GTF文件
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/gencode.vM23.annotation.gtf.gz

基于比對的hisat2流程

構(gòu)建索引
ref=$wkd/reference/mm10.fasta
hisat2-build -p 40 $ref hisat.mm10
# 歷時(shí)：00:14:39

如果不想自己構(gòu)建索引，可以直接下載hisat官網(wǎng)數(shù)據(jù)

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
比對過程
# 以其中SRR7815790為例
index=$wkd/hisat/hisat.mm10
hisat2 -p 10 -x $index -1 SRR7815790_1.fastq.gz -2 SRR7815790_2.fastq.gz  --new-summary -S SRR7815790_hisat.sam

samtools sort -O bam -@ 10  -o SRR7815790_hisat.bam SRR7815790_hisat.sam 

samtools index -@ 10 -b SRR7815790_hisat.bam

比對結(jié)果

# 看到有一個(gè)細(xì)胞（SRR7816020）質(zhì)量不好遮糖，后面需要去除
SRR7815790.hisat.log:   Overall alignment rate: 90.14%
SRR7815850.hisat.log:   Overall alignment rate: 87.49%
SRR7815870.hisat.log:   Overall alignment rate: 87.39%
SRR7815890.hisat.log:   Overall alignment rate: 90.08%
SRR7815910.hisat.log:   Overall alignment rate: 90.32%
SRR7815960.hisat.log:   Overall alignment rate: 91.73%
SRR7815980.hisat.log:   Overall alignment rate: 86.01%
SRR7816020.hisat.log:   Overall alignment rate: 61.82%
SRR7816100.hisat.log:   Overall alignment rate: 88.14%
SRR7816120.hisat.log:   Overall alignment rate: 90.21%
SRR7816130.hisat.log:   Overall alignment rate: 87.30%
SRR7816140.hisat.log:   Overall alignment rate: 89.59%
SRR7816160.hisat.log:   Overall alignment rate: 84.33%
SRR7815790.hisat.log:Execution time for SRR7815790 hisat and sam2bam : 115.350829 seconds
SRR7815850.hisat.log:Execution time for SRR7815850 hisat and sam2bam : 149.416581 seconds
SRR7815870.hisat.log:Execution time for SRR7815870 hisat and sam2bam : 218.976666 seconds
SRR7815890.hisat.log:Execution time for SRR7815890 hisat and sam2bam : 168.738353 seconds
SRR7815910.hisat.log:Execution time for SRR7815910 hisat and sam2bam : 214.020996 seconds
SRR7815960.hisat.log:Execution time for SRR7815960 hisat and sam2bam : 169.597353 seconds
SRR7815980.hisat.log:Execution time for SRR7815980 hisat and sam2bam : 146.707250 seconds
SRR7816020.hisat.log:Execution time for SRR7816020 hisat and sam2bam : 3.484966 seconds
SRR7816100.hisat.log:Execution time for SRR7816100 hisat and sam2bam : 131.300152 seconds
SRR7816120.hisat.log:Execution time for SRR7816120 hisat and sam2bam : 133.382700 seconds
SRR7816130.hisat.log:Execution time for SRR7816130 hisat and sam2bam : 124.014219 seconds
SRR7816140.hisat.log:Execution time for SRR7816140 hisat and sam2bam : 166.373111 seconds
SRR7816160.hisat.log:Execution time for SRR7816160 hisat and sam2bam : 209.587995 seconds
定量
gtf=$wkd/reference/gencode.vM23.annotation.gtf
featureCounts  -T 10 -p -t exon -g gene_name -a $gtf -o $wkd/hisat/count/hisat2_counts.txt $wkd/hisat/*.bam 1>featureCounts.log 2>&1

不基于比對的salmon流程

salmon需要對轉(zhuǎn)錄本構(gòu)建索引绣的。因此只能使用參考轉(zhuǎn)錄組，而不能使用基因組

構(gòu)建索引
ref=$wkd/reference/gencode.vM23.transcripts.fa
salmon index -t $ref -i salmon.mm10
定量
# 以其中SRR7815790為例
index=$wkd/salmon/salmon.mm10
salmon quant -i $index -l A --validateMappings \
        -1 SRR7815790_1.fastq.gz -2 SRR7815790_2.fastq.gz \
        -p 10 -o SRR7815790_quant 1>SRR7815790.salmon.log 2>&1

運(yùn)行速度的確很快

SRR7815790.salmon.log:Execution time for SRR7815790 salmon : 36.016640 seconds
SRR7815850.salmon.log:Execution time for SRR7815850 salmon : 54.478129 seconds
SRR7815870.salmon.log:Execution time for SRR7815870 salmon : 76.499805 seconds
SRR7815890.salmon.log:Execution time for SRR7815890 salmon : 58.166234 seconds
SRR7815910.salmon.log:Execution time for SRR7815910 salmon : 56.277339 seconds
SRR7815960.salmon.log:Execution time for SRR7815960 salmon : 52.918936 seconds
SRR7815980.salmon.log:Execution time for SRR7815980 salmon : 47.993982 seconds
SRR7816020.salmon.log:Execution time for SRR7816020 salmon : 8.981974 seconds
SRR7816100.salmon.log:Execution time for SRR7816100 salmon : 40.505802 seconds
SRR7816120.salmon.log:Execution time for SRR7816120 salmon : 52.856553 seconds
SRR7816130.salmon.log:Execution time for SRR7816130 salmon : 38.405255 seconds
SRR7816140.salmon.log:Execution time for SRR7816140 salmon : 53.823648 seconds
SRR7816160.salmon.log:Execution time for SRR7816160 salmon : 75.848338 seconds

Salmon的結(jié)果和Hisat的不同欲账，它的每個(gè)數(shù)據(jù)都是獨(dú)立的文件夾屡江，其中quant.sf 就是每個(gè)樣本的定量結(jié)果

對于salmon的定量結(jié)果，需要用到R里面的tximport 導(dǎo)入

tximport函數(shù)主要需要兩個(gè)參數(shù)：定量文件files和轉(zhuǎn)錄本與基因名的對應(yīng)文件tx2gene

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

####################
# 配置files路徑
####################
dir <- file.path(getwd(),'quant/')
dir
files <- list.files(pattern = '*sf',dir,recursive = T)
files <- file.path(dir,files)
all(file.exists(files))

####################
# 配置tx2gene
####################
# https://support.bioconductor.org/p/101156/
# BiocManager::install("EnsDb.Mmusculus.v79")
# 如何構(gòu)建赛不？
if(F){
  library(EnsDb.Mmusculus.v79)
  txdf <- transcripts(EnsDb.Mmusculus.v79, return.type="DataFrame")
  mm10_tx2gene <- as.data.frame(txdf[,c("tx_id", "gene_id")])
  head(mm10_tx2gene)
  write.csv(mm10_tx2gene,file = 'mm10_tx2gene.csv')
}
tx2_gene_file <- 'mm10_tx2gene.csv'
tx2gene <- read.csv(tx2_gene_file,row.names = 1)
head(tx2gene)

####################
# 開始整合
####################
library(tximport)
library(readr)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,ignoreTxVersion = T)
names(txi)
head(txi$length)
head(txi$counts)
# 發(fā)現(xiàn)目前counts的列名還沒有指定

####################
# 添加列名
####################
# 獲取導(dǎo)入文件名稱的ID惩嘉，如SRR7815870
library(stringr)
# 先得到SRRxxxx_quant這一部分
n1 <- sapply(strsplit(files,'\\/'), function(x)x[12])
# 再得到SRRxxxx這一部分
n2 <- sapply(strsplit(n1,'_'), function(x)x[1])

colnames(txi$counts) <- n2
head(txi$counts)
             
             
####################
# 操作新得到的表達(dá)矩陣
####################
salmon_expr <- txi$counts
# 表達(dá)量取整
salmon_expr <- apply(salmon_expr, 2, as.integer)
head(salmon_expr)
# 添加行名
rownames(salmon_expr) <- rownames(txi$counts)
dim(salmon_expr)

save(salmon_expr,file = 'salmon-aggr.Rdata')

好，現(xiàn)在有了salmon和featureCounts兩種方法得到的矩陣踢故，就可以來比較一下

rm(list = ls())
options(stringsAsFactors = F)
load('salmon-aggr.Rdata')
dim(salmon_expr)
colnames(salmon_expr)
salmon_expr[1:3,1:3]

# 為了方便比較文黎，我們選第一個(gè)樣本SRR7815790
salmon_SRR7815790 <- salmon_expr[,1]
names(salmon_SRR7815790) <- rownames(salmon_expr)

###########################
# 然后找featureCounts結(jié)果
###########################
ft <- read.table('../hisat/hisat2_counts.txt',row.names = 1,comment.char = '#',header = T)
ft[1:4,1:4]
colnames(ft)
# 兩種取子集的方法，grep返回結(jié)果的序號殿较；grepl返回邏輯值
ft_SRR7815790 <-  ft[,grep(colnames(salmon_expr)[1], colnames(ft))]
ft_SRR7815790 <-  ft[,grepl(colnames(salmon_expr)[1], colnames(ft))]

names(ft_SRR7815790) <- rownames(ft)

###########################
# 取salmon和featureCounts公共基因
###########################
# salmon使用的Ensembl基因耸峭，而featureCounts得到的是Symbol基因
# 先進(jìn)行基因名轉(zhuǎn)換
library(org.Mm.eg.db)
columns(org.Mm.eg.db)
gene_tr <- clusterProfiler::bitr(names(salmon_SRR7815790),
                      "ENSEMBL","SYMBOL",
                      org.Mm.eg.db)
# 要找salmon在gene_tr中的對應(yīng)位置，然后取gene_tr的第二列symbol斜脂，因此match中把salmon放第一個(gè)位置
names(salmon_SRR7815790) <- gene_tr[match(names(salmon_SRR7815790),gene_tr$ENSEMBL),2]

# 找共有基因
uni_gene <- intersect(names(salmon_SRR7815790),names(ft_SRR7815790))
length(uni_gene)

plot(log1p(salmon_SRR7815790[uni_gene]),
     log1p(ft_SRR7815790[uni_gene]))
結(jié)論就是：有些基因用不同的流程檢測效果是不一樣的抓艳，使用不同的參考數(shù)據(jù)得到的結(jié)果也是有區(qū)別的，但整體上一致
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市玷或，隨后出現(xiàn)的幾起案子儡首，更是在濱河造成了極大的恐慌，老刑警劉巖偏友，帶你破解...
    沈念sama閱讀 206,126評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件蔬胯，死亡現(xiàn)場離奇詭異，居然都是意外死亡位他，警方通過查閱死者的電腦和手機(jī)氛濒，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來鹅髓，“玉大人舞竿，你說我怎么就攤上這事×耄” “怎么了骗奖？”我有些...
    開封第一講書人閱讀 152,445評論 0贊 341
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵，是天一觀的道長醒串。 經(jīng)常有香客問我执桌，道長，這世上最難降的妖魔是什么芜赌？ 我笑而不...
    開封第一講書人閱讀 55,185評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任仰挣，我火速辦了婚禮，結(jié)果婚禮上缠沈，老公的妹妹穿的比我還像新娘膘壶。我一直安慰自己，他們只是感情好洲愤，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,178評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布香椎。 她就那樣靜靜地躺著，像睡著了一般禽篱。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上馍惹，一...
    開封第一講書人閱讀 48,970評論 1贊 284
  • 城市分裂傳說
    那天躺率，我揣著相機(jī)與錄音，去河邊找鬼万矾。 笑死悼吱，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的灌旧。 我是一名探鬼主播碟案，決...
    沈念sama閱讀 38,276評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼触创，長吁一口氣：“原來是場噩夢啊……” “哼锯蛀！你這毒婦竟也來了遇西？” 一聲冷哼從身側(cè)響起馅精，我...
    開封第一講書人閱讀 36,927評論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤，失蹤者是張志新（化名）和其女友劉穎粱檀，沒想到半個(gè)月后洲敢，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
    沈念sama閱讀 43,400評論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡茄蚯，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,883評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年压彭，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片渗常。...
    茶點(diǎn)故事閱讀 37,997評論 1贊 333
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡壮不，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出皱碘，到底是詐尸還是另有隱情询一，我是刑警寧澤，帶...
    沈念sama閱讀 33,646評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布尸执，位于F島的核電站家凯，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏如失。R本人自食惡果不足惜绊诲，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,213評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望褪贵。 院中可真熱鬧掂之，春花似錦、人聲如沸脆丁。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評論 0贊 19
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽槽卫。三九已至跟压，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間歼培，已是汗流浹背震蒋。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留躲庄，地道東北人查剖。 一個(gè)月前我還...
    沈念sama閱讀 45,423評論 2贊 352
  • 代替公主和親
    正文 我出身青樓，卻偏偏與公主長得像噪窘，于是被迫代替她去往敵國和親笋庄。 傳聞我的和親對象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,722評論 2贊 345