這次要練習的文章數(shù)據(jù)出自：
Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing. Nat Commun 2018 Dec 4;9(1):5150.
（一）原始數(shù)據(jù)下載
首先，在GEO網(wǎng)站搜索：GSE111229，可以看到這里有768個原始數(shù)據(jù)。
看到圖里的SRA數(shù)據(jù)的SRP133642，copy下來，然后去https://www.ebi.ac.uk/這里網(wǎng)站，這里可以下載SRA文件和fastq文件。
把剛才查到的SRP133642輸入進去，點擊搜索。
得到以下的頁面：
點擊SRP133642（source前面的那個小號的，不是上面的大號字）
你就可以看到700多個原始文件的列表，分了好多頁。這時可以點擊Select columns選擇你想要下載的一些選項，你可以只下載它們的編號，也可以選擇下載它們的下載地址的鏈接，還有你可以選擇下載SRA文件還是fastq文件。這里我只需要它們的fastq下載地址的鏈接，選好之后點擊768 results后面的TEXT，然后下載下來的是一個txt文件，打開后是這樣的：
這是所有文件的fastq文件類型下載的鏈接地址，因為我之后要用ascp下載。
輸入代碼，可以先看一下前幾行：
$ head fastqid.txt
ftp.sra.ebi.ac.uk/vol1/fastq/SRR679/001/SRR6790711/SRR6790711.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR679/002/SRR6790712/SRR6790712.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR679/003/SRR6790713/SRR6790713.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR679/004/SRR6790714/SRR6790714.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR679/005/SRR6790715/SRR6790715.fastq.gz
ftp.sra.ebi.ac.uk/vol1/fastq/SRR679/006/SRR6790716/SRR6790716.fastq.gz
接下來把每一行的開頭部分“ftp.sra.ebi.ac.uk”去掉，這里就需要用到sed。把每一行的開頭部分都替換成空的字符，并且替換好的輸入到一個新的文件里：
$ sed -e 's/ftp.sra.ebi.ac.uk//g' fastqid.txt >> fastqid_trim.txt
再查看一下新的文件的前幾行，是不是開頭部分都去掉了：
$ head fastqid_trim.txt
/vol1/fastq/SRR679/001/SRR6790711/SRR6790711.fastq.gz
/vol1/fastq/SRR679/002/SRR6790712/SRR6790712.fastq.gz
/vol1/fastq/SRR679/003/SRR6790713/SRR6790713.fastq.gz
/vol1/fastq/SRR679/004/SRR6790714/SRR6790714.fastq.gz
/vol1/fastq/SRR679/005/SRR6790715/SRR6790715.fastq.gz
然后就是用ascp批量下載文件了：
$ ascp -QT -l 300m -P33001 -k1 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list fastqid_trim.txt /media/yanfang/FYWD/RNA_seq/fastq_files/single_cell_data/
一些參數(shù)的解釋：
-T： 取消加密，否則有時候數(shù)據(jù)下載不了。
-Q：一般都加，我也不知道啥意思。
-l 設置最大傳輸速度，一般200m到500m。
-P33001： UDP端口。
-k： 斷點續(xù)傳，一般設置為值1。
-i： 提供私鑰文件的地址。就是你的*.openssh這個文件放在哪里了。
--mode：模式，是上傳還是下載。可以選send, recv。
--host：服務器
--user：用戶
--file-list：你要下載的鏈接地址存在了哪個文件里，這個文件的名字。
最后要寫明你下載的這些文件想存放在哪個文件夾里。
NOTE：這里提一句，如果運行ascp的時候出現(xiàn)了Session Stop的報錯，可以嘗試以下解決方法（參考：從NCBI-SRA和EBI-ENA數(shù)據(jù)庫下載數(shù)據(jù))：
第一、把ascp 的-l 不要設置的太大了，一般設置為500m就可以了；（實際上我只用了300m,我這個破電腦。。。）
第二、在root權(quán)限下設置udp端口：
iptables -I INPUT -p udp --dport 33001 -j ACCEPT
iptables -I OUTPUT -p udp --dport 33001 -j ACCEPT
一開始我這兩條都試了，但還是出錯，整的我要崩潰，弄了一個多小時，就是一直報錯，結(jié)果后來發(fā)現(xiàn)，是我最后的輸出文件的途徑寫錯了。。。文件夾的名字少了一個字母。。。好吧，我這個智商真是非常感人了。
寫好命令行，就可以愉快的下載了，這篇文章是單端測序，所以每一個樣品的Fastq文件只有一個：
$ ascp -QT -l 300m -P33001 -k1 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list fastqid_trim.txt /media/yanfang/FYWD/RNA_seq/fastq_files/single_cell_data/
SRR6790711.fastq.gz 100% 12MB 19.9Mb/s 00:06
SRR6790712.fastq.gz 100% 12MB 15.4Mb/s 00:13
SRR6790713.fastq.gz 100% 6453KB 14.1Mb/s 00:18
SRR6790714.fastq.gz 100% 6633KB 18.3Mb/s 00:21
SRR6790715.fastq.gz 100% 8949KB 15.7Mb/s 00:26
SRR6790716.fastq.gz 100% 6274KB 17.1Mb/s 00:31
SRR6790717.fastq.gz 100% 13MB 19.7Mb/s 00:38
SRR6790718.fastq.gz 100% 5129KB 22.5Mb/s 00:41
SRR6790719.fastq.gz 100% 10MB 25.4Mb/s 00:46
...中間的部分我就不貼上來了，700多行呢，最后下載好會有一個顯示：
Completed: 8002176K bytes transferred in 4628 seconds
(14163K bits/sec), in 768 files.
（二）質(zhì)控
理論上要看每個樣品的質(zhì)量如何，決定是否要進行trim，因為這篇文章是跟著“單細胞天地”微信公眾號跑一遍，主要是先大體的了解一下單細胞測序分析的過程，所以這次就沒有做trim，而且我挑了一些樣品做fastqc，發(fā)現(xiàn)總體質(zhì)量還湊合，基本上都是這樣的，而且沒有接頭污染，接頭的含量顯示為0。
（三）比對
在單細胞天地的公眾號文章里（單細胞轉(zhuǎn)錄組上游分析之shell回顧 ），作者用的是hisat2進行比對的。而這篇文獻的作者在方法的部分對測序分析的過程是這樣描述的：“Single-end 43?bp long reads were aligned at ESCG to mm10 mouse genome using STAR v2.3.043 with default settings”?？磥碜髡哂玫氖荢TAR來進行比對的。那么我就順便又查了一下hisat2和STAR有什么區(qū)別。根據(jù)大神的說法（http://www.itdecent.cn/p/479c7b576e6f）：
HISAT2找到junction正確率最高，但是在總數(shù)上卻比TopHat和STAR少。HISAT2的二類錯誤(納偽）比較少，但是一類錯誤（棄真）就高起來。就唯一比對而言，STAR是三者最佳的，主要是因為它不會像TopHat和HISAT2一樣在PE比對不上的情況還強行把SE也比對到基因組上。而且在處理較長的read和較短read的不同情況，STAR的穩(wěn)定性也是最佳的。就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍。
嗯。。?？吹阶詈笠痪?，我還是決定用hisat2吧。。。我這破電腦如果用STAR去比對700多個樣品的話，估計會冒煙。
另外，這篇文獻的單細胞測序是使用了spike-in作為control的。目的是用來消除一些因素對差異分析的影響。文獻中提到：“sequenced using the Smart-seq2 protocol with exogenous RNA controls from External RNA Controls Consortium (ERCC) spiked into the cell lysates”，可以看到作者是使用的spike-in RNA的。如果存在spike-inRNA，這些序列就要在比對前（構(gòu)建基因組index）的時候加上去。在之后的read count步驟里，還需要把spike-in和參考基因組GTF文件合并在一起。因為外源加入的spike-in的物種不是小鼠的，如果不在比對的時候加入到參考基因組里，那么你比對出來所有的spike-in是不會有任何Read比對到參考基因組上的，那么spike-in的意義就沒有了。我下載的是mm10基因組的fasta文件和GTF文件，ERCC(spike-in)的下載地址是：https://www.thermofisher.com/search/results?query=ERCC92.gtf%20&focusarea=Search%20All。合并的代碼(https://assets.thermofisher.com/TFS-Assets/LSG/manuals/cms_095048.txt)：
$ cat ERCC92.fa >> mm10_plus_spikein.fa
$ cat ERCC92.gtf >> gencode.vM22.annotation_plus_ERCC.gtf
合并之后需要自己構(gòu)建參考基因組索引，因為我之前下載的小鼠的索引是沒有spike-in的，現(xiàn)在加上了spike-in，所以要自己構(gòu)建了。代碼很簡單：
$ hisat2-build -p 20 mm10_plus_spikein.fa genome
構(gòu)建索引的速度取決于你的電腦，我這個破電腦用了2個多小時。也還算可以接受的。
下面就開始比對了：
ls *.gz | while read i;do hisat2 -p 10 -x /media/yanfang/FYWD/RNA_seq/ref_genome/index/mm10_plus_ERCC_hisat2/genome -U $i -S /media/yanfang/FYWD/RNA_seq/sam_files/single_cell_sam/${i%%.*}.sam;done
402682 reads; of these:
402682 (100.00%) were unpaired; of these:
55003 (13.66%) aligned 0 times
347439 (86.28%) aligned exactly 1 time
240 (0.06%) aligned >1 times
86.34% overall alignment rate
395135 reads; of these:
395135 (100.00%) were unpaired; of these:
50638 (12.82%) aligned 0 times
344300 (87.13%) aligned exactly 1 time
197 (0.05%) aligned >1 times
87.18% overall alignment rate
212554 reads; of these:
212554 (100.00%) were unpaired; of these:
55914 (26.31%) aligned 0 times
156545 (73.65%) aligned exactly 1 time
95 (0.04%) aligned >1 times
73.69% overall alignment rate
......
參數(shù)解釋：
這里用了一個while語句，關于while語句，可以參考我之前貼上來的文章（結(jié)構(gòu)化命令 之 while命令）。
-P 10：用10個線程。
-x:后面跟的是小鼠基因組index的存放位置。
-U：單端數(shù)據(jù)文件。
-S:輸出文件存放的位置。
%%.*：取輸入原始測序文件(i)的前綴。
（四）sam轉(zhuǎn)換成bam文件
$ ls *.sam | while read i ;do (samtools sort -O bam -@ 10 -o /media/yanfang/FYWD/RNA_seq/bam_files/single_cell_bam/$(basename ${i} ".sam").bam ${i});done
其中一個小tip就是：basename這里，會返回每個$i，也就是每個sam文件的名稱，然后".sam"就是將名稱中的.sam去掉，于是只留下了前綴名，這樣才能更改文件名為bam。（這個方法來自：單細胞轉(zhuǎn)錄組上游分析之shell回顧 ）
（五）構(gòu)建索引
接下來就是構(gòu)建每一個bam文件的索引了～
$ ls *.bam |xargs -i samtools index {}
這里我查了一下xargs是什么命令：
xargs 可以將管道或標準輸入（stdin）數(shù)據(jù)轉(zhuǎn)換成命令行參數(shù)，也能夠從文件的輸出中讀取數(shù)據(jù)。xargs 是一個強有力的命令，它能夠捕獲一個命令的輸出，然后傳遞給另外一個命令。之所以能用到這個命令，關鍵是由于很多命令不支持|管道來傳遞參數(shù)，所以就有了 xargs 命令。一般是和管道一起使用。
命令格式：
somecommand |xargs -item command
（六）count計數(shù)
之前我自己練習的bulk-RNA-seq分析時，這一步read計數(shù)我都是使用的htseq-count這個軟件。而這次看單細胞天地文章的時候，作者在進行read count時使用了featureCounts。原諒我的無知。。。之前沒有聽說過這個軟件，于是上網(wǎng)搜索了featureCounts和htseq-count的區(qū)別：
（1）首先，這兩個軟件運行后read count數(shù)是不一樣的。（例如：htseqcount和featurecount計數(shù)不一樣 ）
（2） HTSeq-count對于多重比對的reads則是采取舍棄策略。（參考轉(zhuǎn)錄組定量可以用替換featureCounts代替HTSeq-count）
（3） Htseq-count慢！如果樣本數(shù)是幾百個時，其所消耗的時間是以天記的。（怪不得我之前6個8個的樣本就要過夜跑。。。也不排除我這破電腦的原因）
（4）不僅支持基因或轉(zhuǎn)錄本的定量，還可以進行exons, promoter, gene bodies, genomic bins and chromosomal locations的定量。
“單細胞天地”給出的代碼是：
# 還是先指定gtf的位置
$ gtf=/YOUR_GTF_PATH/
$ featureCounts -T 5 -t exon -g gene_id -a $gtf -o all.id.txt *.bam 1>counts.id.log 2>&1 &
# 如果是雙端測序，可以加上-p選項
-T：使用的線程數(shù)。
-t：Specify feature type in GTF annotation. `exon' by default. Features used for read counting will be extracted from annotation using the provided value.要計數(shù)的feature名稱，也就是GTF的第三列信息，默認是exon。
-a：Name of an annotation file.指定使用的注釋文件。
-o: Name of the output file including read counts. 輸出文件。
-g ：Specify attribute type in GTF annotation. "gene_id" by default. Meta-features used for read counting will be extracted from annotation using the provided value.提取的GTF最后一列attribute信息，默認是gene_id。如果你不想讓你的count 結(jié)果是ENSMUSG。。。一堆數(shù)字，而是基因的具體名稱，則需要把gene_id換成gene_name。
下面根據(jù)我的文件位置，更改一下代碼：
$ featureCounts -T 5 -t exon -g gene_id -a /media/yanfang/FYWD/RNA_seq/ref_genome/gencode.vM22.annotation_plus_ERCC.gtf -o /media/yanfang/FYWD/RNA_seq/matrix/singlecell/all.id.txt *.bam 1>counts.id.log 2>&1
這里說一句，1>counts.id.log的意思是標準輸出；2>&1 的意思就是將標準錯誤重定向到標準輸出。可以把counts.id.log看作"黑洞"，它等價于一個只寫文件，所有寫入它的內(nèi)容都會永遠丟失，而嘗試從它那兒讀取內(nèi)容則什么也讀不到。
最后得到兩個文件：一個是矩陣，另外一個是矩陣的總結(jié)信息
（七）過濾表達矩陣
得到矩陣以后，先讀取矩陣
> a=read.table('all.id.txt',header = T ,sep = '\t')
單細胞測序嘛，肯定很多基因的表達量都是0，比如說矩陣的第一行：
> x=a[1,]#把第一行賦值給x
> table(x>1)#看一眼第一行里這個基因表達量大于1的細胞有多少個
FALSE TRUE
768 5
#結(jié)果在這么多的細胞里，只有5個細胞有這個基因的表達
那我們有700多個細胞，總不能一個一個的看吧？所以要用到循環(huán)。循環(huán)的具體代碼完全按照這篇文章的來：獲取Github代碼包以及準備工作。這篇文章的代碼也是解析的非常詳細。
>floor(ncol(a)/50)#用總列數(shù)除以50然后向下取整，結(jié)果就是15
也就是說，只要一行中至少要在15個樣本中有表達量。這個對矩陣的過濾，沒有統(tǒng)一的標準。我查詢了幾篇文章，幾篇文章里的過濾標準都不一樣。有的是按照每一個細胞里的基因數(shù)過濾的（有的人要求一個細胞里要有7000個基因表達，有的人100就夠了），有的是按照每一個基因在一定的樣品數(shù)里有表達，這個就不一定了，要看你的樣品量和你的分析目的。我特意貼了三篇對矩陣過濾的文章，可以看出這一步完全看心情？？？比如下面這篇：
使用limma、Glimma和edgeR，RNA-seq數(shù)據(jù)分析易如反掌
“在任何樣本中都沒有足夠多的序列片段的基因應該從下游分析中過濾掉。這樣做的原因有好幾個。 從生物學的角度來看，在任何條件下的表達水平都不具有生物學意義的基因都不值得關注，因此最好忽略。 從統(tǒng)計學的角度來看，去除低表達計數(shù)基因使數(shù)據(jù)中的均值 - 方差關系可以得到更精確的估計，并且還減少了在觀察差異表達的下游分析中需要進行的統(tǒng)計檢驗的數(shù)量?！?/p>
還有這篇：
單細胞RNA-seq分析
1.去除在任何細胞中都不表達的基因。
2.具有少量reads/molecules，很可能是已經(jīng)被破壞或捕獲細胞失敗，應該移除這類細胞
3.手動過濾細胞：大多數(shù)細胞檢測到的基因在7000-10000之間，這對于高深度scRNAseq來說是正常的。但是這也取決于實驗協(xié)議和測序深度。比如說基于droplet或其他測序深度較低的方法，通常其每個細胞檢測到的基因較少。如果細胞檢出率相等，則分布應近似正常，因此我們移除分布尾部的那些細胞（檢測到的基因少于7000的細胞）
4.自動過濾細胞：自動異常檢測來識別可能存在問題的細胞：scater包創(chuàng)建一個矩陣（行代表cell，列代表不同的QC度量），然后PAC提供了按照QC度量排序的cells的2D表示，然后使用來自mvoutlier包的方法檢測異常值。在這里插入圖片描述
5.過濾基因：移除表達水平被認為是undetectable的基因。We define a gene as detectable if at least two cells contain more than 1 transcript from the gene. If we were considering read counts rather than UMI counts a reasonable threshold is to require at least five reads in at least two cells. 然而很多時候閾值取決于測序深度，而且很重要的一點是基因過濾一定要在細胞過濾之后，因為一些基因可能僅僅被檢測到只存在在低質(zhì)量的細胞里。
還有一篇是這樣寫的：
單細胞測序（scRNA-seq）通關||數(shù)據(jù)處理必知必會
根據(jù)基因的表達量等特征，對細胞進行過濾，通常的做法就是指定一個閾值，比如要求一個細胞中檢測到的基因數(shù)必須大于100，才可以進入到下游分析，如果小于這個數(shù)字，就過濾掉該細胞。需要強調(diào)的是，在設定過濾的閾值時，需要人為判斷，這樣的設定方式會受到主觀因素的干擾，所以往往都會指定一個非常小的過濾范圍，保證只過濾掉極少數(shù)的離群值點。
綜上所述，矩陣過濾沒有統(tǒng)一標準，請自行嘗試。
繼續(xù)說矩陣的過濾代碼：
上面知道了 x>1 返回邏輯值0和1，0為FALSE，1為TRUE?，F(xiàn)在我們要找一行中總共有多少TRUE，就用sum計算一下(因為會忽略掉0的影響)
>sum(x>1) > floor(ncol(a)/50)# 當然第一行會返回FALSE，也就表明我們要去掉這一行內(nèi)容
>a[sum(x>1) > floor(ncol(a)/50),]# 就把不符合要求的第一行去掉了
這是對一行的操作，那么對于上萬行的矩陣怎么操作？答：用apply()
> allsample=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]
#首先，對a這個矩陣進行操作
#其次，對行(也就是1表示)進行操作；如果對列操作，用2表示。
#復雜的操作先寫上 function(x){}，這是一個標準格式。
#最后，把篩選完成的矩陣賦值給新的變量。
> dim(allsample)#看看新的矩陣里保留了多少基因？16289個。
[1] 16289 773
(八)過濾后操作
過濾之后，看一下樣品里的spike-in的情況（有關spike-in的介紹，這篇文章里的幾篇文獻有提到單細胞RNA-seq入門文獻學習）：
> grep('ERCC',rownames(allsample))
integer(0)
簡直驚呆了我！明明有92個spike-in的control，怎么一個也沒有？后來我看了一下表達矩陣才發(fā)現(xiàn)：基因名，或者說spike-in的名字，被當成了第一列，并沒有成為行名，所以才檢索不到有spike-in的存在。那么需要把第一列設置成行名：
> rownames(a)<-a[,1] #將數(shù)據(jù)框的第一列作為行名
> b<-a[,-1] #將數(shù)據(jù)框的第一列刪除，只留下剩余的列作為數(shù)據(jù)
然后就可以看spike-in的情況了：
> grep('ERCC',rownames(allsample))
[1] 15928 15929 15930 15931 15932 15933 15934 15935 15936 15937 15938 15939 15940 15941 15942 15943 15944 15945 15946 15947
[21] 15948 15949 15950 15951 15952 15953 15954 15955 15956 15957 15958 15959 15960 15961 15962 15963 15964 15965 15966 15967
[41] 15968 15969 15970 15971 15972 15973 15974 15975 15976 15977 15978 15979 15980 15981 15982 15983 15984 15985 15986 15987
[61] 15988 15989
（九）去除文庫大小差異（參考：常說的表達矩陣，那得到之后呢？）
每個細胞測得數(shù)據(jù)大小不同，這樣是沒辦法看高表達還是低表達的，必須先保證基數(shù)一樣才能比較，cpm（counts per million）這個算法就是做這個事情的。
cpm是歸一化的一種方法，代表每百萬堿基中每個轉(zhuǎn)錄本的count值。
注意：這個算法只是校正文庫差異，而沒有校正基因長度差異。要注意我們分析的目的就是：比較一個基因在不同細胞的表達量差異，而不是考慮一個樣本中不同兩個基因的差異，因為"沒有兩片相同的樹葉”這個差異是正常的。但是同一個基因由于某種條件發(fā)生了改變，背后的生物學意義是更值得探索的。
用起來很簡單，有現(xiàn)成的函數(shù)cpm() ，然后我們再用log將數(shù)據(jù)降個維度，但保持原有數(shù)據(jù)形狀不變：log2(edgeR::cpm(dat)+1) 。
意思就是：cpm需要除以測序總reads數(shù)，而這個值作為分母會導致結(jié)果千差萬別，有的特別大有的很小。為了后面可視化不受極值的影響，用log轉(zhuǎn)換一下可以將數(shù)值變小，并且原來大的數(shù)值最后還是大，并不改變這個現(xiàn)實。
# 先看看前6行6列的數(shù)據(jù)
>finaldata[1:6,1:6]
X.SS2_15_0048_A1 X.SS2_15_0048_A2 X.SS2_15_0048_A3 X.SS2_15_0048_A4 X.SS2_15_0048_A5 X.SS2_15_0048_A6
ENSMUSG00000102269.1 0 0 0 0 0 0
ENSMUSG00000098104.1 5 0 0 0 1 1
ENSMUSG00000103922.1 1 0 1 0 0 0
ENSMUSG00000033845.13 0 0 9 1 19 0
ENSMUSG00000025903.14 202 46 29 80 0 27
ENSMUSG00000033813.15 1 5 1 0 9 8
# 比如先計算一下第五個樣本的總測序量
> sum(dat[,5])
[1] 215937 #說明第五個樣品有這么多的read數(shù)
# 那么對于第5個樣品的第4個基因read數(shù)是19
# 計算它的cpm值：count值*1000000/總測序reads
> 19*1000000/215937
[1] 87.98863
# 和標準公式比較看看，結(jié)果完全相同
> edgeR::cpm(finaldata[,5])[4]
[1] 87.98863
# 因此最后就是
> finaldata2=log2(edgeR::cpm(finaldata)+1)
#把去除文庫大小差異后的matrix保存成一個table
> write.table(finaldata2,"adjust_library_matrix.txt",sep="\t")
