『壹』 RNAseq教程(2.2)
1.Mole 1 - Introction to RNA sequencing
2.Mole 2 - RNA-seq Alignment and Visualization
3.Mole 3 - Expression and Differential Expression
4.Mole 4 - Isoform Discovery and Alternative Expression
5.Mole 5 - De novo transcript reconstruction
6.Mole 6 - Functional Annotation of Transcripts
用HISAT2比對基因組和轉錄組。
首先，為對齊結果創建適當的輸出目錄
HISAT2的輸出是每個數據集的SAM/BAM文件。
參考HISAT2幫助手冊獲得更多說明：
HISAT2基本用法:
額外參數如下:
注意:在上面的對齊中，我們將每個庫視為一個獨立的數據集。如果你有一個庫的多個數據，你可以在一個HISAT2命令中將它們對齊在一起。要組合多個reads文件，您需要為'-1'輸入參數提供所有read1文件作為逗號分隔的列表，然後為'-2'輸入參數提供所有read2文件作為逗號分隔的列表，(其中兩個列表的順序相同):還可以使用samtool
將HISAT2 sam文件轉換為bam文件，並按對齊位置排序
將所有UHR數據和所有HBR數據合並成一個BAM文件。注意:這可以通過幾種方式來完成，比如「samtools merge」，「bamtools merge」，或者使用picard-tools(見下文)。我們選擇第三種方法是因為它在合並bam頭信息方面做得最好。注意:sambamba也保留頭部信息。
計算對齊(BAM)文件，確保所有文件都成功創建(總共應該有8個)
任務:對額外的數據集執行一些比對。用你在上面學到的技巧來對齊閱讀。嘗試使用HISAT2。還要練習將SAM轉換為BAM文件，並合並BAM文件。
在練習3中創建的名為「practice」的單獨工作目錄中進行分析。
.sam和.bam文件之間有什麼區別?
sam文件是一個純文本序列比對映射文件。bam文件是相同信息的二進制壓縮版本。
如果您像上面所做的那樣對結果BAM文件進行排序，那麼結果是否按read名稱排序?還是position?
按照position排序
可以查看BAM文件的哪些列以確定排序的樣式?
第一、第三和第四列包含reads名稱、染色體和位置。
可以使用什麼命令僅查看BAM頭?
『貳』 RNA-seq原理
測序技術發展：
1977Sanger測序--1996焦磷酸測序--2003cmPCR--2003ZMW---2012納米孔測序
RNA-seq的一些技術限制，測序誤差主要由生物學誤差(生物學重復，比如取30隻小鼠采樣)和技術性誤差(技術性重復，比如對1隻小鼠采樣3次)造成，如果想要得到的數據為無偏的，那麼生物學重復最重要，因為生物個體代表著樣本，而技術手段只會造成不可控干擾。總的來說，只做技術性重復的實驗結果偏差最大，技術性重復+生物學重復的實驗結果偏差也可能較大，除非生物學重復遠大於技術性重復(因為當生物學重復次數不足時，技術性重復能擴大樣本單一的影響)，無論如何，多做生物學重復，這有助於你的結論被其他人復現。
原理詳解：
A 為了保證細胞在標記的過程中是單獨分開的，10X開發了微流體設備(microfuidic device)進行預處理，設備有三個上樣孔，分別加入你的1.樣本細胞懸液(Sample) 2.凝膠小球(Beads) 3.分離液(Oil)，下圖為具體設備的示意圖。
當我們把樣本細胞懸液加入設備時，每一個細胞會與凝膠小球單獨結合，然後被分離液包裹，形成一個油包水的密閉小液滴(droplet)。進一步地，細胞和凝膠小球相遇不久後會裂解，釋放出裡面的各種物質，RNA(mRNA、tRNA、rRNA)，蛋白質，脂質，DNA等。實際上Beads上聯接了不同的接頭，其中有一個接頭包含ploy(dT)序列，在細胞裂解後釋放的核酸中，只有mRNA帶有polyA tail，於是Beads的poly(dT)接頭就可以從眾多的裂解產物里捕獲到mRNA(實際上drop-seq採用3'端測序，就是為了檢測polyA tail)。
Master Mix中帶有反轉錄試劑，當mRNA被捕獲後，就可以從它的3『端開始作為模板，進行反轉錄出cDNA的第一條鏈，這第一條鏈就沿著poly(dT)序列延申，長在了beads上，形成了圖一7中的STAMPs，接著我們把反轉錄出來的cDNA序列洗脫，以cDNA的第一條鏈為模板，進行PCR，合成cDNA的第二條鏈，然後就是我們熟悉的cDNA擴增以及illumina測序。
如何確定測序序列來自哪個細胞？single cell的RNA-seq和bulk的RNA-seq的最大區別是什麼？是barcode，或者說是cell barcode(實際上DNA自帶barcode，cell barcode是人為控制的)。每一種single cell的beads上都有著相同的cell barcode(beads與beads間的cell barcode是不同的)，假設每個beads只捕獲一個cell，那麼則每個cell都被cell barcode 單獨標記了。
如何保證每個beads只捕獲一個cell？第一是控制cell和beads的流速，第二是beads的數目遠遠超過cell的數目，即絕大多數的beads都是空的，只有少數的才捕獲到了cell。但是還是有個別的droplet裡面會兩個或者更多的細胞，這就需要進行質控(QualityControl)。
接下來可以參照10X Genomics的說明書詳解single cell RNA-seq的barcode。
實際上beads上一開始只接了Read1、Barcode、Poly(dT)。
名詞解釋：
Poly(dT): 用來和mRNA的polyA結合，捕獲mRNA
UMI: 用來標記不同的PCR產物(用於count計數)。為了減少由於復制引起的誤差(重復抽樣導致重復計數)，人們在一些單細胞測序的步驟中增加了UMI(unique molecular identifiers)，UMIs 是由 4-10 個隨機核苷酸組成的序列，在 mRNA 反轉錄後，進入到文庫中，每一個 mRNA，隨機連上一個 UMI，因此可以計數不同的 UMI，最終計數 mRNA 的數量。
10X Barcode: 用來標記不同的single cell
Sample Index: 用來標記不同的sample
P5和P7: 用來進行illumina的橋式PCR測序
Truseq Read 1、2: 用來進行連接beads，cDNA的PCR擴增和加P7接頭
在這些序列中，P5、P7、Truseq Read 1、2 的序列是已知的。
其他的序列是怎麼一步一步添加上去的？
具體步驟：
利用Poly(dT)來捕獲mRNA，在mRNA的5'端插入TSO(Template Switch Oligo模板切換低聚糖)引物，然後從mRNA的polyA開始反轉錄，直至mRNA的DNA序列被轉錄完成，然後在beads序列的3'端插入CCC，再對mRNA的TSO進行反轉錄，至此完成了cDNA的第一條鏈(序列順序和mRNA逆序)。上述步驟很重要，因為中間cDNA的序列我們是不知道的(儀器測序長度有限)，如果不加上這個接頭，就沒有辦法設計引物來合成cDNA的第二條鏈。
將mRNA溶解，對cDNA的第一條鏈加入UMI引物，以cDNA的第一條鏈為模板合成cDNA的第二條鏈。最後使用PCR(聚合酶鏈式反應)對cDNA(拷貝DNA)進行擴增(為了富集)。
PCR原理
因為II代測序(NGS)的illumina測序不能測很長的seq，約為200-700bp，所以不能測得mRNA全長，因此需要進一步把合成的cDNA利用酶打斷到illumina能測的長度(長度有些隨機，比如300bp的cDNA能通過頭尾150bp完整測序，但700bp的cDNA只能通過頭尾150bp測序+參考基因組推斷出來)。然後在cDNA的3'端插入Truseq Read2引物(和Truseq Read1引物匹配為頭尾，中間序列就是reads)、P5、P7。
最後的測序數據(reads)從Truseq Read1後的10X Barcode開始，一直到Truseq Read2為止。
PCR擴增是對cDNA單鏈進行復制，後面的橋式PCR是對完整的樣本進行復制(增加數據深度)，總的來說各個cDNA呈均勻分布，然後進行抽樣。
RNA-seq plications有PCR plication(最主要)、cluster plication、optical plication。
實際上儀器會對核苷酸進行染色，然後判斷顏色確定ATCG鹼基，因此有很多原因會導致機器誤判，和後續QC有關。
1.某些核苷酸對顏色附著不明顯
2.大片區域顏色相同(相同類型核苷酸)，而其中僅有幾個顏色不同的點(不同類型的核苷酸)
『叄』 單細胞測序（sc-RNA seq）分析：Seurat4.0系列教程、高級分析、文章復現
時代的洪流奔涌而至，單細胞技術也從舊時王謝堂前燕，飛入尋常百姓家。雪崩的時候，沒有一片雪花是無辜的，你我也從素不相識，到被一起捲入單細胞天地。R語言和Seurat已以勢如破竹之勢進入4.0時代，天問一號探測器已進入火星烏托邦平原了，你還不會單細胞嗎？那麼為了不被時代拋棄的太遠，跟著我們一起通過學習seurat系列教程入門單細胞吧。
首先我們學習經典的標准流程，這個教程我當初入門時候曾經花1000購買過，也曾花6000購買過，不同的單細胞培訓班，買的是一樣的教程。現在免費送給你，別說話，開始學吧！
Seurat4.0系列教程1：標准流程 - (jianshu.com)
Seurat4.0系列教程2：多模式數據聯合分析 - (jianshu.com)
Seurat4.0系列教程3：合並兩個樣品的10x數據集 - (jianshu.com)
Seurat4.0系列教程4：整合分析 - (jianshu.com)
Seurat4.0系列教程5：交互技巧 - (jianshu.com)
Seurat4.0系列教程6：常用命令 - (jianshu.com)
Seurat4.0系列教程7：數據可視化方法 - (jianshu.com)
Seurat4.0系列教程8：細胞周期評分和回歸分析 - (jianshu.com)
Seurat4.0系列教程9：差異表達檢測 - (jianshu.com)
Seurat4.0系列教程10：降維 - (jianshu.com)
Seurat4.0系列教程11：使用sctransform - (jianshu.com)
Seurat4.0系列教程12：大數據集整合的方法 - (jianshu.com)
Seurat4.0系列教程13：使用互惠 PCA （RPCA） 快速整合數據 - (jianshu.com)
Seurat4.0系列教程14：整合scRNA-seq and scATAC-seq數據 - (jianshu.com)
Seurat4.0系列教程15：映射和注釋查詢數據集 - (jianshu.com)
Seurat4.0系列教程16：多模式參考映射注釋細胞 - (jianshu.com)
Seurat4.0系列教程17：Mixscape - (jianshu.com)
Seurat4.0系列教程18：Weighted Nearest Neighbor Analysis（WNN)） - (jianshu.com)
Seurat4.0系列教程19：多線程並行策略 - (jianshu.com)
Seurat4.0系列教程20：單細胞對象的格式轉換 - (jianshu.com)
Seurat4.0系列教程21： 結合Cell Hashing分析雙細胞 - (jianshu.com)
Seurat4.0系列教程22：空間轉錄組的分析、可視化與整合 - (jianshu.com)
Seurat4.0系列教程告一段落，但這決不是終點。這個系列教程只是給大家打開一扇窗，知道Seurat4.0有這些功能可用，少走彎路，起到一個拋轉引玉的作用，後續還要自己深入研究。不要像我當初入門單細胞之時，為了找到整合方法耗費大量時間及不必要的金錢。君子生非異也，善假於物也。學，然後知不足。加油吧！少年！用好單細胞測序及分析這個技術，為人類健康研究做出自己的貢獻！不足之處，懇請批評指正！------以此紀念粉絲突破1000。
CellChat三部曲1：使用CellChat對單個數據集進行細胞間通訊分析 - (jianshu.com)
CellChat 三部曲2：使用CellChat 對多個數據集細胞通訊進行比較分析 - (jianshu.com)
CellChat 三部曲3：具有不同細胞類型成分的多個數據集的細胞通訊比較分析 - (jianshu.com)
RNAvelocity系列教程1：RNA速率簡介及scVelo安裝 - (jianshu.com)
RNAvelocity系列教程2：使用 Seurat 和 scVelo 估計 RNA 速率 - (jianshu.com)
RNAvelocity系列教程3：使用Seurat和velocyto估算RNA速率 - (jianshu.com)
RNAvelocity系列教程4：velocyto.R的使用教程 - (jianshu.com)
RNAvelocity系列教程5：scVelo 的基本知識 - (jianshu.com)
RNAvelocity系列教程6：# scVelo用於RNA 速率基本流程 - (jianshu.com)
RNAvelocity系列教程7：# scVelo實戰-胰腺數據集 - (jianshu.com)
RNAvelocity系列教程8：# scVelo實戰-齒狀回數據集 - (jianshu.com)
RNAvelocity系列教程9：# scVelo應用-動力學模型 - (jianshu.com)
RNAvelocity系列教程10：# scVelo應用-微分動力學 - (jianshu.com)
使用fgsea進行單細胞轉錄組GSEA富集分析 - (jianshu.com)
Seurat24節氣之21大雪---vars.to.regress = "percent.mt"排除線粒體基因 - (jianshu.com)
Seurat24節氣之2雨水---resolution參數 - (jianshu.com)
單細胞轉錄組基礎分析五：細胞再聚類 - (jianshu.com)
單細胞測序文章圖表復現01—讀入文件的表達矩陣構建Seurat對象 - (jianshu.com)
單細胞測序文章圖表復現02—Seurat標准流程之聚類分群 - (jianshu.com)
單細胞測序文章圖表復現03—區分免疫細胞和腫瘤細胞 - (jianshu.com)
單細胞測序文章圖表復現04—聚類分群的resolution參數問題 - (jianshu.com)
CNS圖表復現11—RNA-seq數據可不只是表達量矩陣結果 - (jianshu.com)
『肆』 單細胞綜述之整合分析
文章發表於nature review genetics： Integrative single- cell analysis ，作者是Tim Stuart與 Rahul Satija 。做過單細胞分析的對他們應該不陌生。
scRNA-seq技術的發展契合了研究個體細胞表觀遺傳、空間研究、蛋白質組與譜系信息的方法需要，這為研究多類型數據的綜合方法提出了獨特的機遇與挑戰。綜合分析可以發現細胞之間的模式關系，獲取細胞的整體狀態信息，產生涵蓋不同樣本與不同研究手段的數據集。該文重點討論了單細胞基因表達數據與其他類型的單細胞分析方法的整合。
多模態（Multimodal）數據 ：多種類型數據的組合，如RNA與蛋白質數據組合，是一種多維度數據，類似多組學。
單模態 ：單個類型數據
Pseudotime ：擬時分析
聯合聚類（Joint-clustering） ：通過聯合不同類型數據對細胞進行分組。
典型相關分析（CCA） ： 利用綜合變數對之間的相關關系來反映兩組指標之間的整體相關性的多元統計分析方法。
動態時間規整（Dynamic time warping） ：一種局部拉伸或壓縮兩個一維矢量以校正一個矢量相對於另一個矢量的滯後的方法。
MNNs ：標准化基因表達空間中最臨近的細胞。聚類用校正批次效應。
梯度推進（Gradient boosting） ：一種預測模型演算法。
隨著分子生物學、微流控與納米技術的發展，催生了許多類型的單細胞測序技術。過去的方法集中在單模態測量上，如DNA序列、RNA表達量和 染色質可及性 上。雖然這些技術促進了我們對細胞多樣性與發育景觀的理解，但是它們並不能很好地解析單細胞內分子間互作關系。而這些互作關系是深入探索細胞狀態的關鍵。此外，隨著可用數據集規模的快速增長，迫切需要用於標准化與聯合分析且考量到批次效應與個體差異的計算方法。
scRNA-seq是應用最為廣泛的單細胞測序技術之一。而後出現了一系列互補技術如單細胞基因組、表觀基因組和蛋白質組分析技術，涵蓋了單細胞基因組測序（ Vitak, S. A. et al., 2017 ; Navin, N. et al., 2011 ）、染色質可及性（ Pott, S., 2017 ; Corces, M. R. et al., 2016 ; Buenrostro, J. D. et al., 2015 ; Cusanovich, D. A. et al., 2015 ; Lake, B. B. et al., 2018 ）、DNA甲基化（ Luo, C. et al., 2017 ; Smallwood, S. A. et al., 2014 ; Guo, H. et al., 2013 ; Mulqueen, R. M. et al., 2018 ）、膜蛋白（ Stoeckius, M. et al., 2017 ; Peterson, V. M. et al., 2017 ）、小RNA（ Faridani, O. R. et al., 2016 ）、組蛋白修飾（ Gomez, D. te al., 2013 ; Rotem, A. et al., 2015 ）和染色體構象（ Ramani, V. et al., 2017 ; Nagano, T. et al., 2013 ）等技術。目前已開發出研究單細胞空間結構和譜系信息的方法（ Frieda, K. L. et al., 2017 ; Shah, S. et al., 2016 ）。
單細胞多模態綜合分析方法示意
單模態與多模態分析方法匯總
CEL-seq ：線性擴增測序法
CITE- seq ：膜蛋白豐度與基因表達水平測定
G＆T-seq ：基因組轉錄組測序
LINNAEUS ：譜系追蹤
MARS-seq ：大規模平行單細胞RNA測序
MEMOIR ：譜系與空間結構測定
MERFISH ：主要是細胞間結構測定
osmFISH ：環狀單分子熒光原位雜交，空間結構測定
REAP- seq ：膜蛋白豐度與基因表達水平測定
scATAC-seq ：單細胞空間結構測定
scBS-seq ：單細胞甲基化測序
scChIP-seq ：單細胞ChIP-seq
scGESTALT ：結合CRISPR-cas9的譜系追蹤弄方法
scHi-C-seq ：測定染色體組裝
sciATAC-seq ：結合index轉座酶的scATAC-seq
sci-CAR ：利用index聯合分析mRNA和染色質可及性譜
sci-MET ：利用index分析單細胞甲基化水平
sci-RNA-seq ：結合index的scRNA-seq
SCI-seq ：單細胞組合標記測序，檢測CNV
scM＆T-seq ：單細胞甲基化組和轉錄組測序，可研究未知的DNA甲基化與基因表達之間的關系
scNOMe- seq ：核小體佔位與甲基化組測序
scRRBS ：單細胞限制性代表區域甲基化測序
scTHS- seq ：單細胞轉座體超敏性位點測序
seqFISH ：內含子序貫熒光原位雜交，擴展觀測到基因數量
snmC-seq ：單核甲基胞嘧啶測序
SNS ：單核測序
SPLiT-seq ：丐版scRNA-seq
STARmap ：原位單細胞測序
理想的實驗流程應當全面洞悉細胞的所有方面，包括分子狀態、空間構象、胞外環境互作的全部過程。盡管當下技術手段無法做到，但多模態技術與綜合計算方法可以是我們離該目標越來越近。文章希望提出整合單細胞轉錄組學、基因組學、表觀組學與蛋白組學的數據統一分析方法，重點在結合其他數據類型分析scRNA-seq數據，尤其是整合來自於同一細胞的不同類型數據。
文章分為四大塊，首先探討了多模態單細胞分析方法，其次研究了不同實驗不同數據整合分析，然後討論了單細胞空間測序數據整合分析方法，最後給出了整合分析方法的前景與必要性。
最初的單細胞分析方法主要關注細胞某狀態下的某類分子水平。而現在更引人矚目的是同時分析單細胞內多種分子以建立更全面的單細胞分子視圖。通常這些方法是將scRNA-seq數據與其它分析手段的結合，目前主要有四種策略從單細胞中得到多模態數據：
嚴格來說這種方法算單模態。
一些scRNA-seq workflow採用流式分選細胞，隨後進行scRNA-seq（MARS-seq/Smart-seq/2），這樣可以同時獲得單細胞與對應的熒光信號，將熒光所表示的蛋白質水平與轉錄組在同一細胞中關聯（ Ramsköld, D. et al., 2012 ; Jaitin, D. A. et al., 2014 ; Picelli, S. et al., 2013 ）。早期研究（ Hayashi, T. et al., 2010 ）利用FACS結合半定量RT-PCR（作者稱之為FBSC‐PCR），結合scRNA-seq，明確了細胞表面marker可以區分細胞類型與狀態( Wilson, N. K. et al., 2015 ;該文結合了Smart-seq2)，( Paul, F. et al., 2015 ;該文結合了MARS-seq)和鑒定稀有細胞的思路。 Paul, F. et al., 2015 與 Nestorowa, S. et al., 2016 利用該workflow研究發現了小鼠造血祖細胞由轉錄組定義不同細胞簇的免疫表型， Wilson, N. K. et al., 2015 則分離了小鼠HSCs，鑒定細胞維持乾性相關的表面marker。但是囿於熒光光譜的重疊現象，利用該法測到的每個細胞的參數范圍有限。
針對熒光無法分選的部分，FACS顯然是不合適的，尤其是需要同時測得單細胞基因組與胞內蛋白的scRNA-seq實驗。此時需要物理分離或通過不同tag篩選出不同組分。
G&T-seq通過加入oligo(dT)特異性分離mRNA同時保留基因組DNA從而實現了基因組轉錄組平行測序（ Macaulay, I. C. et al., 2015 ）DR-seq通過則通過加入barcode特異擴增cDNA序列實現基因組轉錄組平行測序（ Dey, S. S. et al., 2015 ）。這使得單細胞基因表達水平與其對應基因型聯系起來，深度揭示單細胞間DNA拷貝數變異與染色體重排對下游mRNA豐度的具體關聯。這些方法適用於研究體細胞基因高度變異的腫瘤組織。
DNA甲基化與轉錄組水平結合研究是基於 Macaulay, I. C. et al., 2015 的G&T-seq和 Smallwood, S. A. et al., 2014 的scBS- seq技術發展的，同普通BSP一樣，用亞硫酸氫鈉處理DNA片段隨後進行擴增，結合G&T-seq，可以分析同一細胞內的DNA甲基化模式和基因表達數據（ Angermueller, C. et al., 2016 ）。由於DNA甲基化存在不穩定性和異質性，因此若要研究DNA甲基化與基因表達間的關系，則必須將表觀基因組變異與細胞間的異質性區別開來。
通過DNA甲基化與轉錄組關聯分析，為啟動子甲基化與基因表達間的負相關性提供深層次的證據。此外，利用barcode系統選擇性標記基因組DNA與cDNA，結合index系統，可以對數千個單細胞進行染色質可及性與基因表達水平間的關聯分析，同時鑒定出影響基因表達的順式調控元件（ Cao, J. et al., 2018 ）。
關於胞內蛋白與mRNA關聯研究，有兩種思路可供借鑒。其一（ Darmanis, S. et al., 2016 ）是將FACS sort到的細胞裂解後分離裂解液，分別進行蛋白質與RNA定量。作者採用 PEA (鄰近探針延伸分析) 檢測蛋白並用RT-qPCR定量，採用qRT-PCR定量mRNA。該法可以同時檢測82個mRNA/75個蛋白；其二（ Genshaft, A. S. et al. ）是將FACS sort到的細胞在微流控晶元中同時進行逆轉錄和PEA而不分離裂解液。該法可以同時檢測96個mRNA/38個蛋白。這兩種方法檢測的蛋白與mRNA數量與質量均有限。
這些技術的出現表明若將可以細胞信息轉化為有序的barcode，我們就可以在分析單細胞轉錄組時將這些信息同時獲取。這種策略不僅適用於分析細胞的自然狀態，也適用於大規模基因擾動研究。目前有Perturb-Seq（ Dixit, A. et al., 2016 ）和CRISPR-Seq（ Adamson, B. et al., 2016 ; Datlinger, P. et al., 2017 ; Jaitin, D. A. et al., 2016 ），他們將scRNA-seq與CRISPR-cas9結合進行遺傳篩選，使得研究正向遺傳學的大規模基因擾動試驗成為可能。具體原理是給單個基因擾動和受到影響的細胞添加barcode，通過scRNA-seq能夠鑒定出這兩者，從而推斷CRISPR靶向基因和由此產生的單個細胞的轉錄譜間的關系。目前應用在基因調控網路（ Dixit, A. et al., 2016 ）、未折疊蛋白反應（ Adamson, B. et al., 2016 ）、免疫細胞分化發育（ Datlinger, P. et al., 2017 ）和T細胞受體激活（ Jaitin, D. A. et al., 2016 ），非編碼區調控元件（ Klann, T. S. et al., 2017 ）。此外，還可以結合CRISPR-dcas9系統，擴展到轉錄調控、表觀遺傳調控領域中（ Thakore, P. I. et al., 2016 ; Liu, X. S. et al., 2016 ; Hilton, I. B. et al., 2015 ; Konermann, S. et al., 2015 ; Gilbert, L. A. et al., 2017 ）,18年發展了同時靶向和敲除基因的技術（ Boettcher, M. et al., 2018 ）。
另一個應用是結合CRISPR-cas9的譜系追蹤技術。單細胞譜系追蹤是去年的大熱方向之一，此處提到三種mRNA+lineage方法： scGESTALT 、 ScarTrace 、 LINNAEUS 。這三種方法各有不同，但大體是利用CRISPR-cas9連續切割結合到基因組上的barcode，細胞會用NHEJ來應對這種損傷。但NHEJ容易出錯，從而在DNA序列中產生隨機突變，這些突變通過細胞分裂進行遺傳，結合scRNAseq利用這些突變作為復合barcode來構建組織或器官發育譜系。
另一種略有不同的方法是 MEMOIR ，它結合smFISH與CRISPR-cas9系統，可以同時檢測細胞譜系與空間位置。
普通的scRNA-seq流程除了可以做轉錄本豐度外，還可以進行諸如體細胞突變、遺傳變異、RNA isoform等分析。
關於體細胞突變目前已有研究（ Lodato, M. A. et al., 2015 ）,該文通過對人大腦的少量單細胞全基因組測序，分析了發生的細胞突變，構建了人大腦神經細胞譜系。作者發現突變大多發生在高轉錄活性相關位置，這表明可能可以通過scRNA-seq數據來分析神經細胞突變情況，根據轉錄狀態重構神經細胞譜系。此外，分析scRNA-seq數據中的拷貝數變異，可以研究癌症非整倍體與異質性等情況（ Tirosh, I. et al., 2016 ; Fan, J. et al., 2018 ）。
單細胞分析也為理解DNA自然變異如何影響基因表達與細胞狀態提供了新思路。有研究結合GWAS+scRNAseq，鑒定出了不同個體之間的eQTL（ Kang, H. M. et al., 2018 ）。
多模態測序策略正在催生與之相匹配的數據分析方法。多模數據集可以檢測到細胞間的細微差異，而單模數據很可能無法做到這一點。由於scRNAseq數據存在dropout，故而它更容易忽略細胞間的細微差別；但與來自同一細胞的其他數據互補分析可以改善這一問題。例如，很難通過scRNA-seq數據區分不同的T細胞亞群，但聯合膜蛋白分析則可以顯著提高亞群解析度（ Stoeckius, M. et al., 2017 ），同樣，RNA+chromatin、RNA+methylation聯合可能揭示單個細胞間的調控異質性，不再贅述。
單細胞多模態分析思路很可能受到bulk-seq多組學聯合分析的啟發（ Meng, C. et al., 2016 ）， Argelaguet 開發了一種名為MOFA( multi- omics factor analysis)的方法，該方法在多組學bulk-seq數據中效果良好，同時測試了單細胞DNA甲基化數據與RNA數據聯合處理情況，效果也可以。這暗示適用於bulk-seq的多組學數據處理方式可能也適用於單細胞多模態數據。鑒於單細胞數據規模遠超bulk-seq，多視圖機器學習不失為一種重要的補充手段（ Colomé- Tatché, M. & Theis, F. J., 2018 ）。
單細胞多模態研究策略為解析細胞內不同組分間的關系提供了新方法。如CITE-seq和REAP-seq可以輕易鑒別出相關度較低的RNA-protein模塊，表明此處存在活躍的轉錄後調節。還有一個很有意思的是通過測量剪接過的成熟RNA與未剪接RNA的相對豐度，可以建立RNA與蛋白的關聯動態模型（ La Manno, G. et al., 2018 ）。
此外，還可以在不同類型數據間建立統計模型。前面提到的sci-CAR文章建立了染色質可及性與基因表達水平間的統計模型，通過染色質可及性數據估計細胞內基因表達水平（ Cao, J. et al., 2018 ），另一組研究人員建立了gRNA與基因表達水平間的線性回歸模型，用以識別細胞應答的前後關系，重構轉錄網路（Perturb-Seq（ Dixit, A. et al., 2016 ））。通過這種手段可以研究目標物種復雜的調控網路。
前面主要講了在同一測序實驗同一批細胞進行的多模態數據整合，而不同測序實驗數據整合分析才是亟需解決的關鍵問題。同bulk seq 數據一樣，處理批次效應是綜合分析不同實驗室、不同workflow產出數據的首要問題（SVA包（ Leek, J. T. 2014 ））。然而目前bulk seq水平的處理方法無法處理單細胞數據（（ Haghverdi, L, et al., 2018 ，作者用MNN處理數據，該法在 mnnpy 中得到改進）; Butler, A, et al,. 2018 ）。目前最新方法利用 CCA / MNN 可以識別出兩個數據集間共有的部分，判定細胞間共有的生物學狀態,然後以這些相同狀態的細胞為基準消除批次效應。
此處作者介紹了他自己在Seurat V2中開發的方法（ Satija, R, et al., 2015 ;），該法用 CCA 鑒別出不同數據集間相同的細胞類型且可以避免出現由批次效應或常規PCA造成的假陽性細胞類型；接下來採用動態時間規整演算法校正數據集間細胞密度差異。這兩步驟可以將細胞投影到一個低維空間，具有相同生物學狀態的細胞相互接近且消除了不同數據集帶來的影響。
另一種方法即mnnCorrect，最早用於計算機領域圖形識別。該法尋找不同數據集間最接近的細胞，將之判定為潛在的狀態相同細胞，隨後利用成對MNNs距離計算一個批次參數（batch vector），用以校正原始表達矩陣（ Haghverdi, L., 2018 ）。
CCA/mnnCorrect在整合處理不同來源的scRNA-seq數據時表現良好。這將極大提升發現稀有細胞、微弱轉錄差異細胞及與之對應maker的能力（ Haghverdi, L, et al,.2018 ； Butler, A,et al,. 2018 ） 。這為建立一個統一的單細胞參考數據集提供了依據。在此基礎上，scRNA-seq數據整合分析得到了快速發展（ Hie, B. L, et al., 2018 ; Barkas, N. et al., 2018 ; Park, J.-E., 2018 ; Korsunsky, I. et al., 2018 ; Stuart, T. et al., 2018 ; Welch, J. et al., 2018 ）。這種多數據集整合分析的應用遠不止用於校正批次效應這么單一。它可以在單細胞尺度上深入比較細胞間的狀態，發現細胞對環境及基因擾動的特異性響應，對不同疾病及不同治療下的患者的測序數據進行標准化。
scRNA-seq數據整合分析還可以擴展至跨物種分析。 Karaiskos,N 比較了兩種果蠅早期胚胎的空間基因表達模式，通過構建空間基因表達圖譜，該研究系統比較了兩個果蠅的同源基因表達譜，鑒定出了彼此間的進化波動。 Tosches 比較了爬行動物與哺乳動物腦細胞間的相關性。 Baron 分析了人與小鼠胰島細胞scRNA-seq數據，鑒定出了二者間的保守亞群。 Alpert 開發出了cellAlign，在一維水平上比對了人與小鼠的擬時軌跡，發現人胚胎合子激活要比小鼠晚，小鼠中比人活躍的基因皆與蛋白合成相關。跨物種分析未來是光明的，但對於多物種整合分析而言，精確鑒定物種間同源基因是多物種整合分析至關重要的一步。
以細胞分類信息的形式串聯不同的scRNA-seq數據集，或者借鑒到自己實驗中，是優於合並數據集然後de novo聚類這種方法的。且隨著 有參細胞圖譜 的開發，這種方式將更加尋常。目前已開發對應方法： scmap- cell & scmap- cluster ，其中scmap-cell 用乘積量化（ proct quantization ）演算法進行比對，而scmap-cluster則用於識別未知數據集中的cluster。
利用已有的注釋數據集，目前開發的新方法採用 奇異值分解 、 線性判別分析 和 支持向量機 演算法來對細胞進行分類。此外，隨著引用數據集的大小、范圍與深度越來越高，監督聚類在解析細胞類型方面要比無監督聚類強得多。通過以上這些方法，可以更精確地識別並解析細胞亞群。
satija已有相關文章研究： Comprehensive Integration of Single-Cell Data
這一部分講的是將scRNA-seq數據與其它不同來源和類型數據諸如甲基化、染色質結構等整合分析的方法。
將scRNA-seq數據與其它類型、不同來源的單細胞數據整合分析是無法提取到數據間的共同特徵的，因為它們不是一個類型的數據，需要不同的分析方法。這點在基於基因組的數據（如染色質可及性與甲基化數據）與基於基因的數據（如基因與蛋白表達數據）間整合分析尤為明顯。但如果這些數據來自於同一類細胞群，由於存在著共同的生物學狀態，此時可以聯立分析以發現不同數據集類型間的對應關系。
MATCHER 是一種在一維水平上比較不同類型測序數據擬時軌跡的方法。簡單來說就是比對不同類型測序數據的擬時軌跡，以確定這些數據集間的對應關系。這種方法可以識別不同數據集間的「等效細胞」而不需預先知道彼此間的對應關系。開發者用scM&T- seq（ Angermueller, C. et al., 2016 ）和scRNA-seq數據做了驗證，准確預測了DNA甲基化與基因表達之間的關系。
其他sc-seq數據不同於scRNA-seq數據一樣可以藉助Marker解析細胞類型，因此可以利用scRNA-seq解析出的細胞信息為其他sc-seq數據分析做參考。有研究（ Lake, B. B. et al., 2018 ）對不同腦組織切片進行了單核RNAseq（snRNA-seq）與單細胞轉座子超敏性位點測序（scTHS-seq），通過梯度推進演算法利用單細胞基因表達譜指導了染色質可及性測序數據集的細胞分類：作者首先鑒別出snRNA-seq數據集與scTHS-seq數據集共有的細胞亞群，訓練一個可以將基因表達與染色質可及性數據關聯的模型；然後利用該模型去分類scTHS-seq中剩餘未被分類的細胞。這種方法可以更細致地對大腦組織中的細胞進行分類。同樣，可以整合scATAC-seq數據集來分析單細胞DNA甲基化或轉座酶染色質可及性間的細胞分類。
目前正在開發的新方法有利用假定等價特徵、或識別在所有類型數據中的假定相關共享特徵來進行數據交叉模態分類。 Welch 開發了一種集成非負矩陣分解（iNMF）的方法，名為LIGER，可以跨模態整合數據。他們對同一類型 皮質細胞 分別進行了亞硫酸鹽測序（snmC- seq）與scRNA-seq並對其進行了分類。他們假設基因體甲基化與其表達水平負相關從而整合了不同模態測序數據進行細胞分類。在seurat v3.0中，作者也引入了假定等價特徵或關聯特徵進行多模態整合數據細胞分類的方法。這些方法優點如上所述，即可以利用scRNA-seq的細胞分類信息來指導scATAC-seq數據細胞分類，鑒別出染色質可及性與DNA甲基化的細胞特異模塊。
組織中細胞的空間結構常反映出細胞間的功能差異與細胞命運和譜系的差異。不同基因表達引導細胞向不同方向分化，不同細胞精確排列形成不同組織。關鍵是單細胞實驗通常在分析前細胞已被解離，組織原位信息無法保留，scRNA-seq得到的表達譜不能完全反應細胞空間信息。具有相似基因表達譜的細胞可能存在於不同的空間位置中，故而細胞分離過程中空間信息的缺失是很多單細胞實驗的主要缺點。結合高解析度基因表達譜與空間表達圖譜 (spatial expression maps) 將細胞空間坐標與基因表達譜聯系起來，可以解決這一問題。有兩類方法：計算模型或者RNA原位定量，可以同時收集到細胞空間坐標與基因表達值。
『伍』 RNA-seq 分析之我見（一）
先說下生物體內RNA的大致組成：
編碼RNA:根據中心法則我們知道，DNA轉錄為mRNA，mRNA通過tRNA翻譯為蛋白質，蛋白質行使生命功能,例如呼吸，運動，消化等等。人類只有2萬左右個蛋白質編碼基因，這些編碼基因只佔人類全基因組的2%左右。mRNA占細胞RNA總量的2%～5%， tRNA占細胞RNA總量的15%左右。
非編碼RNA:有些DNA轉錄為RNA後，不繼續編碼蛋白質，這種RNA叫非編碼RNA（ncRNA），包括microRNA,lncRNA,cirRNA，之前人們認為這些RNA是「垃圾」，但是近年來的研究證明，這些RNA對編碼基因發揮著重要的調控作用，是當下研究的熱點。
rRNA：核糖體RNA,佔RNA總量的80%左右。
廣義上說占總RNA95%左右的rRNA和tRNA也屬於非編碼RNA，但是一般研究中，使用的是它的狹義概念，即除去rRNA和tRNA之外的非編碼RNA。
正常情況下，非編碼RNA調控基因的轉錄翻譯，這些都是有序進行的。
但是當處於異常條件下，或者由於自身衰老變異或者受到外部的刺激，比如細菌病毒的感染，射線照射等，這之後往往導致非編碼RNA表達的變化，進而影響蛋白表達的變化，從而引起一系列的病理反應，最終導致疾病。
那麼反過來，如果我們想了解某一疾病具體的發病機理，我們是不是可以提取某一疾病狀態下組織或者細胞的總RNA，去分析它們和正常組表達的異同，我們有理由相信，這些差異表達的RNA分子，很可能與發病機制有關，研究這些差異分子，可以給我們對這一疾病的發病機制的研究提供重要線索，從而研發出更有效的診斷和治療方法。
通過上面的分析，接下來面臨的問題就是，我怎麼分析某一疾病狀態下組織或者細胞所有RNA的表達情況，一個一個分析，肯定不現實，而且可能還有很多未被發現但是很重要的分子。怎麼辦？只有一個辦法，轉錄組測序，即RNA-Seq, 某一條件下所有轉錄出來的RNA鹼基序列，我都給你測出來是什麼。
那麼這涉及6個步驟
1、提取組織或細胞總RNA後，除去佔大部分的rRNA和tRNA，剩下編碼RNA 和非編碼RNA
2、對這些RNA進行測序，理想情況下，是直接檢測，但是不現實，只有通過鹼基互補配對的合成過程，才能知道原來樣品中模板的序列，但是這個合成的長度是有限制的，所以只能先把這些RNA切割成小片段，再檢測這些小片段的序列。具體原理見陳巍學基因視頻。這個過程得到兩種數據，一種是許許多多的鹼基序列，一個是這些序列的表達頻率。也就是一個是RNA是什麼鹼基序列，一個是RNA表達了多少量
3、由於上一步把RNA切割了，好像是一塊拼圖打散了，所以，這一步需要將這一個個的小塊再重新拼成一個完整的圖片。也就是比對，將檢測到的RNA鹼基序列，比對到參考基因組上，看某段RNA位於參考基因組的哪段序列上。這一步就好像一個拼好的拼圖，上面有高高低低的小塊，有些分子表達量高，它對應那個小塊就高，反之就低。通過這一步，實驗組和對照組都得到一個高高低低的拼圖。
4、把實驗組和對照組的拼圖比較一下，看哪些RNA小塊表達量是不一樣的。或者你高我低，或者我高你低，從而得到這些差異表達基因名字的列表。因此這一步的結果都是一些基因名字或者轉錄本編號了。
5、將這些差異表達的分子，進行下游功能分析，比方看看它們都跟什麼信號通路相關啊，可能跟什麼功能有聯系啊等等。這一步得到的就是很多結果圖了。
6、下一步就是將篩選到的差異基因，結合你感興趣的生物學功能或者過程，挑選出幾個，進行再進一步的機制研究。這步就是濕實驗了，也是決定文章層次的核心，這是需要人力和財力，再加上運氣的事情，不過就算不做這一步，前5步也能發篇小文章灌灌水了。由於這一步涉及基礎醫學的機制研究方法，不在這篇文章的討論范圍內。
樣品送測序儀器後，也就是上述第2步後得到會產生大量的數據，可能是多少個G的級別，有幾萬，甚至幾十萬的鹼基序列，首先你要比對到參考基因組吧，然後你要看看實驗組和對照組哪些基因表達有差異吧，其實這不是很復雜的事情，就是數據量太大了，如果就幾十個，你完全可以用EXCEL查找，再標記，但是幾十萬個基因，誰能做到啊。所以現在需要一種工具，可以對數據進行批量編輯和操作。
感謝計算機發達的技術，前人早就幫我們想出來。Linux操作系統就可以實現對大量數據的批量編輯
，R語言可以實現大量數據的統計和做圖。
好了，我們的下一步就是學習Linux操作系統和R語言了。
但是這兩個部分包含了很多很多的知識，我們完全零基礎，要是從頭開始學，效率有點低，畢竟不是專業計算機出身，不需要一下子把所有東西都學會，先把目前需要掌握的學到，將來再舉一反三，慢慢學其它的。
所以現在就開始模擬實戰，從一個測序數據的樣本開始，看看是怎麼一步步得到文章中的結果的。
未完待續...
這兩天寶寶得了幼兒急疹，耽誤了幾天，現在真的是上有老，下有小了，生活的壓力會逼得你迅速成長起來。到這個時候才能深刻體會到時間是非常寶貴的，盡量少干不必要的事情，抓緊提升吧~
『陸』 bulk RNA-Seq （4）合並表達矩陣
上一步定量，我們是對每一個樣本進行的，現在需要將它合並到一個矩陣。這里通過一個腳本merge_to_matrix.pl 來合並。
輸入：每個樣本的定量結果
輸出：
reads count 矩陣： genes.counts.matrix #原始矩陣
表達量矩陣：genes.TMM.EXPR.matrix #標准化後的矩陣
『柒』 WGCNA（轉載）
WGCNA原理及應用
WGCNA介紹：
WGCNA（weighted gene co-expression network analysis，權重基因共表達網路分析）是一種分析多個樣本基因表達模式的分析方法，可將表達模式相似的基因進行聚類，並分析模塊與特定性狀或表型之間的關聯關系，因此在疾病以及其他性狀與基因關聯分析等方面的研究中被廣泛應用。
WGCNA演算法是構建基因共表達網路的常用演算法（詳解： http://www.jianshu.com/p/94b11358b3f3 ）。WGCNA演算法首先假定基因網路服從無尺度分布，並定義基因共表達相關矩陣、基因網路形成的鄰接函數，然後計算不同節點的相異系數，並據此構建分層聚類樹(hierarchical clustering tree)，該聚類樹的不同分支代表不同的基因模塊(mole)，模塊內基因共表達程度高，而分屬不同模塊的基因共表達程度低。最後，探索模塊與特定表型或疾病的關聯關系，最終達到鑒定疾病治療的靶點基因、基因網路的目的。在該方法中mole被定義為一組具有類似表達譜的基因，如果某些基因在一個生理過程或不同組織中總是具有相類似的表達變化，那麼我們有理由認為這些基因在功能上是相關的，可以把他們定義為一個模塊(mole)。這似乎有點類似於進行聚類分析所得到結果，但不同的是，WGCNA的聚類准則具有生物學意義，而非常規的聚類方法(如利用數據間的幾何距離)，因此該方法所得出的結果具有更高的可信度。當基因mole被定義出來後，我們可以利用這些結果做很多進一步的工作，如關聯性狀，代謝通路建模，建立基因互作網路等。
WGCNA的用處：
這類處於調控網路中心的基因稱為核心基因（hub gene），這類基因通常是轉錄因子等關鍵的調控因子，是值得我們優先深入分析和挖掘的對象。
在網路中，被調控線連接的基因，其表達模式是相似的。那麼它們潛在有相似的功能。所以，在這個網路中，如果線條一端的基因功能是已知的，那麼就可以預測線條另一端的功能未知的基因也有相似的功能。
下面的問答來自基迪奧，也能加深對WGCNA的理解
問1、調控網路和共表達網路有什麼區別？
答：調控網路是個更廣泛的概念，而共表達網路是調控網路的一種。
理論上我們可以利用各類信息構建調控網路（表達相關性，序列靶向關系、蛋白互作關系），另外調控網路構建的信息既可以來源真實的實驗驗證的關系，也可以來源生物信息的預測。而共表達網路特指利用基因間的表達相關性預測基因間調控關系的方法，而WGCNA又是共表達網路分析中最有效的方法之一。
問2、WGCNA分析適合的生物物種范圍有規定么？
答：沒有限制。對於任何物種中心法則都是存在的，調控關系對於任何物種都是存在的，所以WGCNA沒有物種限定。
問3、同一物種，不同來源的轉錄組數據（比如不同文章/資料來源的），可以放在一起做WGCNA分析嗎？
答：只要樣本間有相似的生物學意義，是可以合並在一起做分析的。但要注意，不同批次之間的樣本是有批次效應的，所以可能會帶來一些誤差，但是是可以放在一起分析的。
問4、相同材料不同處理之間，可以放在一起做WGCNA分析嗎？比如重金屬和鹽鹼處理。
答：可以的。這也正式WGCNA強大的地方，其可以將不同處理的樣本，合並在一起做分析。其他方法則不一定有這么強大的能力，比如做基因表達趨勢分析時，如果樣本涉及到多個處理不同時期的時候，就不好合並分析（或合並後難以解讀）。但WGCNA的方法關注的是調控關系，所以不管是多少個處理組，都可以很好的整合在一起做分析。
問5、不同批次的數據能放一起做WGCNA嗎？
答：可以的。雖然有批次的干擾，但是干擾對WGCNA網路沒有太大影響。因為WGCNA不是做差異分析，而是基因的共表達。因為批次效應理論上不影響相關性。
問6、不同類型的材料，比如親本和F1，適合放一起進行WGCNA么？
答：如果是一個作圖群體，當然親本與F1是可以放在一起分析的，因為你只關心基因的表達模式，所以把親本加進來是沒有問題的。
問7、沒有生物學重復，共3組，每組5個時間點能夠做嗎？
答：理論上有15個樣本，是可以做WGCNA分析的。並且，分析出來的結果對你的研究應該是非常有用的。至少他會比趨勢分析更有意義，更加准確。
問8、一般說WGCNA的樣品不少於15個，15個樣品考慮重復嗎？不同倍性的材料呢？
答：15個樣本這個是包含了生物學重復，比如5個時間點3個重復；在RNA-seq裡面建議不要用不同倍性材料加進來。除非是有參考的多倍體，如果是無參的多倍體，不同倍性之間差異太大，會讓調控網路不準確。所以用單一倍性的材料做調控網路會更加准確。
問9、可以將RNA-seq數據與蛋白組數據，甲基化數據放一起做WGCNA分析？
答：不能與蛋白數據一起分析。因為WGCNA是基於相關系數的演算法。所以最好一起分析的數據變異度是類似的，RNAseq變異非常大，而蛋白的數據變異很小，兩者的變化不在一個數量級上面。所以兩種數據放在一起分析不合理。
但RNA數據可以嘗試跟甲基化數據一起分析。當然我們也建議RNA數據與代謝組數據一起分析，因為代謝組的數據變異也非常大。
問10、表達量和表達的基因數目差異太大的樣品可以一起分析嗎？比如樣品A有2k個gene表達 而樣品B有2w個gene表達了 AB可以一起分析嗎？
答：做WGCNA分析的時候，不能脫離生物學意義，既然要分析調控網路，那麼應該分析有相似生物學意義的一組基因，比如說拿相似組織來一起做分析，比如不應該拿大腦的樣本與腳趾的樣本合並在一起做分析，因為很顯然，這兩個組織沒有關聯。如果兩個樣本之間是有相關聯的生物學意義，哪怕表達的基因數不一樣，或表達模式差異很大，那依然可以放在一起分析；但如果樣本之間完全沒有生物學意義，那麼分析就沒有意義。
問11、實驗設計是case3個時間點（各點都有三個重復），control同樣的3個時間點（每點三個重復），WGCNA怎麼做？3個時間點和case-control兩個因素能同時考慮進來分析嗎？
答：可以的。做WGCNA是更加合理的，因為有兩個梯度的樣本，如果只是做差異分析的話，邏輯可能非常復雜，做WGCNA分析是對樣本特性更好的解析，可以直觀看到基因在六個處理組裡面是怎樣表達的。
問12、可以拿混合樣本分析嗎？比如一個病原細菌跟人類細胞的基因，能說明細菌跟人類細胞基因有調控關系嗎？
答：可以。前提是病原菌有足夠的數據並定量准確，並且這個分析是非常有意義的，最後可以說明這些病原菌可以調控哪些宿主基因。
問13、但是病原宿主混合分析的話，宿主蛋白不能分泌到宿主體內豈不是WGCNA生物學上也沒有意義嗎？
答：依然有意義。即使病原的基因沒有分泌到宿主裡面，但是病原的蛋白是會影響宿主基因的調控的，比如某個細菌感染某個植物，雖然細菌的蛋白不能直接分泌到植物體內，但會影響植物蛋白的分泌。混在一起分析依然是有意義，可以看到植物裡面到底哪個基因對細菌蛋白產生應答作用。
問14、晶元數據兩分類，每組20個樣本，能否每組單獨做WGCNA？
答：可以。WGCNA還有一種重要功能是做兩個網路的比較，比如病人20個樣本做一個調控網路，健康人做一個調控網路，然後兩個網路做比較。
問15、WGCNA可以用來分析lncRNA對下游基因的調控分析嗎？
答：可以。WGCNA網路有利於預測lncRNA的潛在功能。
問16、構建網路是用所有表達基因還是差異基因？
答：這個是具體問題具體分析。如果使用所有的基因分析，會導致運算量非常大。而也不是所有的基因在這個實驗中都有生物學意義，所以我們會提前做一些過濾。
但用於分析的基因不一定是差異表達基因，有時可以用差異表達基因做一個並集，或通過計算變異系數將變異系數低的基因以及低表達的基因去除。但注意，如果你有關心的特定目標基因的話，應該盡量給予保留。
問17、關注某一個pathway上的基因以及調控因子之間的相關性，構建WGCNA網路的時候屬於這個pathway的基因數量太少會不會影響結果呢？
答：這不是問題。在一個調控網路裡面，樣本的某個pathway上，並不是所有基因參與調控（或存在差異性），所以在做WGCNA分析的時候，會做一些過濾，將有變化的基因挑出來再做分析。即分析的是某個pathway上有變化的基因，不需要分析pathway上所有的基因，只需要分析那些變化的基因就夠了。
問18、前期篩選的時候，要選出在所有樣本中變異系數比較大的基因呢？還是直接用差異表達的基因取並集？用基因還是轉錄本，哪個好呢？
答：兩則都可以，我推薦使用變異系數，選擇那些變異較大的基因，來做下面的分析。然後建議用基因不要用轉錄本，因為轉錄本的定量是不準確的。
問19、變異系數一般取多大？
答：具體問題具體分析。例如，沒有特定目標的時候，可以先計算變異系數，將變異系數的百分之前50來做分析，把變異系數偏低的後面一半過濾掉。
問20、輸入數據用FPKM合適嗎？
答：可以。
問21、RNA seq數據是RSEM值怎麼辦？
答：RSEM值原始輸出結果為reads數，如果是RSEM值建議做一個RPKM校正再做分析。
問22、除了RPKM值以外，做WGANA是否還需要其他數據？TCGA數據可否來做WGCNA分析？
答：在做WGCNA分析必須要用表達量數據，但TCGA的數據某些層級沒有表達量數據，沒有表達量數據自然就無法做WGCNA分析。
問23、請問輸入的基因樣本的矩陣的時候，要不要對數據標准化？
答：做WGCNA分析的時候，不需要對數據進行標准化，輸入RPKM值就足以做這個分析。雖然一些文章會做log2處理，但我認為取了LOG2後，會讓一些表達關系沒有那麼豐富。
問24、每個樣本有3個生物學重復，不需要對三個重復的表達量求平均值代表該樣本嗎？
答：注意，做WGCNA的時候每個樣本是獨立的，三個生物學重復樣本是全部導入做分析，不是取均值再做分析，每個樣本都是獨立的。
問25、如果3個生物學重復，做WGCNA的時候是取三個值，還是用cuffdiff處理後取一個值？
答：如果是生物學重復樣本進行調控網路分析，每個樣本獨立使用，而不是取均值。
問26、請問將樣本信息同模塊特徵值進行相關性分析的時候，樣本信息是怎麼處理的呢？比如不同取樣點、不同性別什麼的，這不是數量性狀信息的，這種情況應該怎麼處理呢？
答：樣本的任何信息都可以做模塊相關性分析。比如相關時間點，可以按照先後量化為12134567。又如不同性別，男與女，可以定義為1，-1。任何性狀量化為數字後，都可以進行相關性分析。
問27、怎麼將模塊與性狀對應起來呢有些性狀不好量化，如果直接將模塊與分組對應，如何實現， 不需要量化指標么？
答：首先需要將性狀量化，如果無法將性狀量化，那麼就無法分析。至於分組信息，也可以量化為類似00001111000（1代表一種組別，2代表另一組組別），實現分組信息的數字化。
問28、基因數量為3w左右時，moles數量為多少結果較為理想？怎麼評價聚類效果的好壞？
答：moles數量沒有標准，moles數量無法評估模塊分的好壞，分組是否合理應該看樹的樹形圖，比如樹的分支很清晰就說明模塊式清晰的。moles數量數由生物性狀決定的。比如樣本表達信息很豐富的時候，moles數量會很多；如果樣本的基因表達相對單一，moles數量就會比較少。
問29、我運行例子的時候，得出來基因之間的direction全是undirected，這和前面的幾種關系有什麼區別？
答：WGCNA是一個undirected的方法，它的網路是無方向的，有相關關系但是無方向。
問30、如果做有向網路的構建，您推薦那些方法？
答：很多方法，例如貝葉斯的方法。
問31、非模式物種可以得出基因之間的相互關系類型么？得出的結果也是undirected么？
答：WGCNA是基於表達兩處理的，所以即使是非模式生物，當然也可以他們之間關系，並且關系也是一個無向網路。
問32、選擇幾個表型數據進行結合分析比較好
答：越多越好，看實驗設計。
問33、感染小鼠，5個時間點，3個重復，找不到合適的表型怎麼辦？
答：如果找不到合適表型，可以找某個時間點應答的基因，本身基因的表達趨勢已經有某種生物學意義的。沒有找到合適表型，也可以看變化趨勢。不一定要做表型的相關分析，其他分析也是很有趣的。例如，可以對模塊功能的富集分析，其實都是可以幫助你找到特定模塊的。所以不用糾結於做某個表型的關聯分析。
問34、weight就是tom值嗎？
答：是的。
問35、剪模塊是怎麼做的？是根據TOM劃分嗎？需要自己設定，還是R自動的？
答：剪模塊是R中自動完成的，不需要劃分，但合並的時候你可以設定一個指標，比如差異度是0.25。
問36、看WGCNA說明是用相異矩陣D（D=1-TOM）去做聚類，然後動態剪切？
答：用TOM值來構建矩陣，TOM值就是兩個樣本的相似度，1-TOM值就是兩個樣本的差異度，相似度與差異度可以理解為一個東西，並不矛盾。
問37、模塊特徵值和樣本性狀相關分析的具體方法是？
答：R包用的是計算相關系數的方法。
問38、WGCNA裡面一般會提到hubgene，如何確定hubgene?
答：在WGCNA分析裡面，每個基因都會計算連通性，連通性高的就是hubgene。
問39、在R中安裝「」WGCNA「」說不適合R3.3.1，那適合哪個版本？
答：WGCNA應該是所有版本都適合，如果版本沒有可以考慮降低R軟體的版本，這個對分析沒有影響。因為不同R版本是一樣的。
問40、用STEM分析的時候擬合多少個模型合適？
答：建議不要超過20個。模塊太多不好分析。
參考網站：
http://tiramisutes.github.io/2016/09/14/WGCNA.html
http://www.jianshu.com/p/94b11358b3f3
http://www.omicshare.com/class/home/index/classdetail?id=20
『捌』 8.單細胞 RNA-seq 聚類分析：整合
目標：
挑戰：
建議：
通常，在決定是否需要執行任何對齊之前，我們總是在 沒有集成的情況下 查看我們的聚類。 不要總是因為您認為可能存在差異而總是執行集成 - 探索數據。 如果我們在 Seurat 對象中一起對兩種條件進行歸一化並可視化細胞之間的相似性，我們就會看到條件特定的聚類：
通常，在決定是否需要執行任何對齊之前，我們需要先查看沒有整合的集群。 不要認為可能存在差異而執行整合——探索數據 。如果我們在Seurat對象中同時對這兩個條件進行歸一化，並可視化細胞之間的相似性，我們就會看到特定條件下的聚類:
不同條件下的細胞聚類，為了確保相同細胞類型的細胞聚集在一起，我們需要在不同條件下整合細胞。
為什麼讓相同細胞類型的細胞聚集在一起很重要？
在所有樣品/條件/模式中我們要鑒別出 細胞類型 ，並在於我們的數據該細胞類型存在，我們需要探索我們的集群之間差異表達的基因。如果我們在一個簇中有我們的對照 T 細胞，而在另一個簇中有我們的受激 T 細胞，那麼當我們試圖在 T 細胞中找到富集的基因時，我們將返回許多作為 T 細胞標記的基因，以及相關基因處理和控制之間的差異。這些可能難以區分，使細胞類型識別更加困難。如果我們在同一簇中擁有相同細胞類型的所有細胞，它還可以通過識別 保守 的細胞類型來幫助我們識別細胞類型基因在簇之間差異表達。此外，如果您想對這些 T 細胞進行任何下游分析（例如控制/治療之間的 DE 分析、配體-受體分析等），重要的是相同細胞類型的細胞存在於相同的簇中.
在本課中，我們將介紹不同條件整合樣本，它改編自 Seurat v3 引導集成教程 。
如果細胞按樣本、條件、批次、數據集、模式進行聚類，則此整合步驟可以大大地改進聚類和下游分析 。
為了整合，我們將使用來自每個組共有的高度可變基因（使用 SCTransform 識別），然後，我們將「整合」或「協調」這些組以覆蓋具有相似「共同生物學特徵」的細胞或組。例如，我們可以整合：
整合是一種強大的方法，它使用這些最大差異的共享資源來跨條件或數據集識別共享的亞群 [ Stuart 和 Bulter 等人。(2018) ]。整合的目標是確保一個條件/數據集的細胞類型與其他條件/數據集的細胞類型一致(例如，控制巨噬細胞與受刺激巨噬細胞一致)。。
具體地說，這種整合方法期望在組內的至少單個細胞子集之間有「對應」或 共享的生物狀態 。集成分析的步驟如下圖所示:
運行的不同分析步驟如下：
現在，使用我們的 SCTransform 對象作為輸入，讓我們執行跨條件的集成。
首先，我們需要指定我們要使用 SCTransform 識別的所有 3000 個最大突變的基因進行整合。默認情況下，此函數僅選擇前 2000 個基因。
現在，我們需要 准備 用於整合的 SCTransform 對象 。
現在，我們將 執行 CCA，找到最好的夥伴或錨點並過濾不正確的錨點 。對於我們的數據集，這最多需要 15 分鍾才能運行。 另外，請注意控制台中的進度條將保持在 0%，但要知道它實際上正在運行。
最後，我們可以 跨條件整合 。
整合後，為了可視化整合的數據，我們可以使用降維技術，例如 PCA 和 Uniform Manifold Approximation and Projection (UMAP)。雖然 PCA 將確定所有 PC，但我們一次只能繪制兩個。相比之下，UMAP 將從任意數量的主要 PC 獲取信息，以在這個多維空間中排列細胞。它將在多維空間中獲取這些距離，並在二維空間中繪制它們，可以保持局部和全局結構。這樣， 細胞之間的距離代表了表達的相似性 。如果您想更詳細地探索 UMAP， 這篇文章 是對 UMAP 理論的一個很好的介紹。
要生成這些可視化結果，我們需要首先運行 PCA 和 UMAP 方法。讓我們從 PCA 開始。
我們可以通過 PCA 映射看到 PCA 很好地覆蓋了這兩個條件。
現在，我們還可以使用 UMAP 進行可視化。運行該方法並繪制圖表。
有時，如果我們 在 conditions 之間拆分繪圖， 更容易查看所有細胞是否對齊，我們可以通過向 DimPlot() 函數添加 split.by 參數實現：
由於整合可能需要一段時間，因此 保存整合 seurat 對象 是一個很不錯的選擇。
『玖』 RNA-Seq歸一化問題
RNA-Seq在衡量基因表達水平時，若單純以比對到基因上的reads數來計算表達量在統計學上是不合理的。影響因素有：
1.基因長度：需要基因長度來比較同一細胞內不同基因之間的表達。RNA-seq實驗中眾所周知的固有技術效果與基因長度有關：RNA（或cDNA）分子在測序之前先進行片段化，較長的轉錄本會比較短的轉錄本被剪切成更多的片段。因此，轉錄本的reads數不僅與其表達水平成正比，而且與其長度成正比。如下圖，我們不能單純的數比對到基因上的read數來比較表達量高低，需要考慮基因長度的影響。
這樣來說，序列長的基因永遠會被認為表達量較高，從而錯誤估計基因真正的表達量。為了消除基因長度產生的固有技術誤差，在過去十年中，已針對RNA-seq數據開發了許多歸一化方法，其中常用的有RPKM、TMM、RLE、upper quartile上四分位處理等。
所以需要對原始的表達矩陣進行標准化，去除掉測序深度和基因長度所帶來的噪音。
是指所有細胞中均要穩定表達的一類基因，其產物是對維持細胞基本生命活動所必需的，管家基因是一類始終保持著低水平的甲基化並且一直處於活性轉錄狀態的基因。
矯正的思路很簡單，就是在變化的樣本中尋找不變的量，那麼在不同RNA-seq樣本中，那些是不變的量呢？一個很容易想到的就是管家基因。但其實這種辦法有一個非常強的先驗假設：管家基因的表達量不怎麼發生變化。其實管家基因有幾千個，這幾千個基因有一定程度上的變化是有可能的。因此這種方法不準確。
在RNA-Seq建庫的過程中摻入一些預先知道序列信息以及序列絕對數量的內參。這樣在進行RNA-Seq測序的時候就可以通過不同樣本之間內參（spike-in）的量來做一條標准曲線，就可以非常准確地對不同樣本之間的表達量進行矯正
數值概念：計算公式：CPM=C/N*1000000
設C為比對到geneA的read 數（read count ）,N為比對到所有基因的總read 數
用途：在某些情況下，只想了解每個基因被覆蓋到的相對read數，而不希望對其做長度校正，就會使用這個指標。CPM 只對read count相對總read 數做了數量的均一化。如果想進行基因間表達量的比較，則不得不考慮基因長度的不同。如果進一步做長度的均一化，就得到了下面的RPKM.
RPKM，全稱為reads per kilobase per million mapped reads，指的是每一百萬個map 上的reads 中，map到外顯子的每1K個鹼基上的read數。
一般來說，基因越長，讀取的次數（深度）越多，自然其有效讀數就越多。而RPKM 就是為了消除這兩個干擾的因素。以更好的比較不同結果。
數值概念：計算公式：RPKM=(1000000 C)/(N L/1000)
設C為比對到geneA的read 數（read count ）,N為比對到所有基因的總read 數，L為gene A的鹼基數，RPKM法能消除基因長度和測序量差異對計算基因表達量的影響，計算得到的基因可直接用於比較不同樣品間的基因表達差異
用途：用於與基因表達量相的後續分析，用於單端測序
計算步驟：
首先對總值數據進行標准化（當然正常來說map到的count 肯定不止這么多，這里只是除了10，但一般而言百萬級的count 可以除一百萬）
相比於RPKM，FPKM 計算的是fragments，也就是一對reads。與RPKM 的差別主要體現在，FPKM在一對reads map上的情況下只計數1，而RPKM 會計為2。適用於雙端測序。
與 RPKM/FPKM 的差別在於，TPM 首先進行了基因長度的標准化，接著再進行了測序深度的標准化。
步驟：1.長度標准化
2.深度標准化
自然從這點來說TPM 的使用范圍更為廣泛。
參考 excel演示DESeq2歸一化原理 - (jianshu.com)
要回答這個問題，我們需要先撇開所有形式上的計算，重新思考到底什麼是RNA轉錄本的表達豐度這個問題。
事實上，對於任何一個制備好測序文庫的待測樣本，它上面任何一個基因的表達量（或者說豐度），都將已是一個客觀存在的值，這個值是不管你改變了多少測序環境都不會變的。而且總共有多少個不同的基因在這這一刻進行了表達，實際上也已經是客觀定好了
此刻，我們可以假定，對於樣本X，其有一個基因g被轉錄了n次，同時樣本X中所有基因的轉錄總次數假定是m次, 那麼正確描述基因g轉錄豐度的值應該是：rg=n/m
同時，樣本X中其他基因的轉錄豐度的計算也應該和上述公式類似。除了要把分子n換為其他基因對應的轉錄次數之外，分母m都一樣。於是有趣的事情就是，所有基因轉錄本豐度的均值mean將是一個恆定不變的數，由以上定義這個數就是：mean = (a+b+c+...+n+...)/m/x
x代表基因個數
此時發現 (a+b+c+...+n+...)=m，公式化簡變為mean = 1/x
這個期望值竟然和測序狀態無關！僅僅由樣本中基因的總數所決定的。也就是說，對於同一個物種，不管它的樣本是哪種組織（正常的或病變的），也不管有多少個不同的樣本，只要它們都擁有相同數量的基因，那麼它們的mean都將是一致的。這是一個在進行比較分析的時候，非常有意義的恆等關系。
由於上面的結果是在理論情況下推導出來的，實際上我們無法直接計算這個r，那麼我們可以嘗試通過其他方法來近似估計r，只要這些近似統計量可以隱式地包含這一恆等關系即可
實際數據來證明FPKM和RPKM犯的錯:
假定有兩個來自同一個個體不同組織的樣本X和Y，這個個體只有5個基因，分別為A、B、C、D和E它們的長度分別如下：
我們可以得到，樣本X和Y的轉錄本的不變數，mean值是:
我們以FPKM的計算的為例子，以下這個表格列出的分別是樣本X和Y在這5個基因分別比對上的fragment數和各自總的fragment數量：
於是，按照以上公式我們可以得到樣本X和Y在這5個基因上的FPKM值分別為：
樣本X在這5個基因上的FPKM均值FPKM_mean = 5,680;
樣本Y在這5個基因上的FPKM均值FPKM_mean = 161,840
首先，我們可以把FPKM的計算式拆分成如下兩部分。
第一部分的統計量是對一個基因轉錄本數量的一個等價描述（雖然嚴格來講也沒那麼等價）:
第二部分的統計量是測序獲得的總有效Fregment數量的百萬分之一:
FPKM其實就只是這兩部分的商！這有道理嗎？分開來看它們似乎都有點道理，但是合起來的時候其實很沒邏輯。
尤其是第二部分（N/10^6），本來式子的第一部分是為了描述一個基因的轉錄本數量，那麼正常來講，第二部分就應該是樣本的轉錄本總數量（或至少是其總數量的等價描述）才能形成合理的比例關系，而且可以看出來FPKM/RPMK是有此意的，這本來就是這個統計量的目的。
可是，它卻失敗了！
N/10^6的大小其實是由RNA-seq測序深度所決定的，並且是一個和總轉錄本數量無直接線性關系的統計量——N與總轉錄本數量之間的關系還受轉錄本的長度分布所決定，而這個分布往往在不同樣本中是有差異的！
比如，在有些基因中，雖然有效比對到基因的Fragment數是相等的，但一般來說長度越長的基因，其被轉錄的次數就越少。那也就是說，N必須將各個被轉錄的基因的長度考慮進去才能正確描述總體的轉錄本數！而FPKM/RPKM顯然沒有做到這一點，這便是FPKM和RPKM出錯的內在原因。
『拾』 RNA-seq比對到轉錄組
如果有注釋較好的參考基因組的物種，可以直接將RNA-SEQ的數據比對到基因組上，也可以比對到從基因組注釋出來的基因上，或者先比到基因組上再預測基因（也就是經典的tophat-cufflinks-cuffdiff流程）。
若沒有注釋較好的參考基因組，建議先用轉錄組序列從頭組裝（推薦trinity），再把從頭組裝的序列作為參考轉錄本，將RNA-SEQ數據比對到從頭組裝的序列中。或者把從頭組裝的序列和對應物種的NCBI unigene的序列整合在一起，去冗餘之後再作為參考轉錄本。