這次要練習的文章數(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權限下設置udp端口：
iptables -I INPUT -p udp --dport 33001 -j ACCEPT
iptables -I OUTPUT -p udp --dport 33001 -j ACCEPT
一開始我這兩條都試了逛裤，但還是出錯，整的我要崩潰猴抹，弄了一個多小時带族，就是一直報錯，結果后來發(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。
（三）比對
在單細胞天地的公眾號文章里（單細胞轉錄組上游分析之shell回顧 ）摹恰，作者用的是hisat2進行比對的辫继。而這篇文獻的作者在方法的部分對測序分析的過程是這樣描述的：“Single-end 43?bp long reads were aligned at ESCG to mm10 mouse genome using STAR v2.3.043 with default settings”∷状龋看來作者用的是STAR來進行比對的姑宽。那么我就順便又查了一下hisat2和STAR有什么區(qū)別。根據(jù)大神的說法（http://www.reibang.com/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晦雨，這些序列就要在比對前（構建基因組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
合并之后需要自己構建參考基因組索引，因為我之前下載的小鼠的索引是沒有spike-in的覆旱，現(xiàn)在加上了spike-in蘸朋，所以要自己構建了。代碼很簡單：
$ hisat2-build -p 20 mm10_plus_spikein.fa genome
構建索引的速度取決于你的電腦通殃，我這個破電腦用了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語句堕担，可以參考我之前貼上來的文章（結構化命令 之 while命令）。
-P 10：用10個線程曲聂。
-x:后面跟的是小鼠基因組index的存放位置霹购。
-U：單端數(shù)據(jù)文件。
-S:輸出文件存放的位置朋腋。
%%.*：取輸入原始測序文件(i)的前綴齐疙。
（四）sam轉換成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轿塔。（這個方法來自：單細胞轉錄組上游分析之shell回顧 ）
（五）構建索引
接下來就是構建每一個bam文件的索引了～
$ ls *.bam |xargs -i samtools index {}
這里我查了一下xargs是什么命令：
xargs 可以將管道或標準輸入（stdin）數(shù)據(jù)轉換成命令行參數(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則是采取舍棄策略昭雌。（參考轉錄組定量可以用替換featureCounts代替HTSeq-count）
（3） Htseq-count慢复唤！如果樣本數(shù)是幾百個時，其所消耗的時間是以天記的烛卧。（怪不得我之前6個8個的樣本就要過夜跑佛纫。。总放。也不排除我這破電腦的原因）
（4）不僅支持基因或轉錄本的定量呈宇，還可以進行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 結果是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)容則什么也讀不到糖耸。
最后得到兩個文件：一個是矩陣，另外一個是矩陣的總結信息
（七）過濾表達矩陣
得到矩陣以后丘薛，先讀取矩陣
> a=read.table('all.id.txt',header = T ,sep = '\t')
單細胞測序嘛嘉竟，肯定很多基因的表達量都是0，比如說矩陣的第一行：
> x=a[1,]#把第一行賦值給x
> table(x>1)#看一眼第一行里這個基因表達量大于1的細胞有多少個
FALSE TRUE
768 5
#結果在這么多的細胞里洋侨，只有5個細胞有這個基因的表達
那我們有700多個細胞舍扰，總不能一個一個的看吧？所以要用到循環(huán)希坚。循環(huán)的具體代碼完全按照這篇文章的來：獲取Github代碼包以及準備工作边苹。這篇文章的代碼也是解析的非常詳細。
>floor(ncol(a)/50)#用總列數(shù)除以50然后向下取整裁僧，結果就是15
也就是說个束，只要一行中至少要在15個樣本中有表達量。這個對矩陣的過濾聊疲，沒有統(tǒng)一的標準茬底。我查詢了幾篇文章，幾篇文章里的過濾標準都不一樣售睹。有的是按照每一個細胞里的基因數(shù)過濾的（有的人要求一個細胞里要有7000個基因表達桩警，有的人100就夠了），有的是按照每一個基因在一定的樣品數(shù)里有表達昌妹，這個就不一定了捶枢，要看你的樣品量和你的分析目的。我特意貼了三篇對矩陣過濾的文章飞崖，可以看出這一步完全看心情烂叔？？固歪？比如下面這篇：
使用limma蒜鸡、Glimma和edgeR胯努，RNA-seq數(shù)據(jù)分析易如反掌
“在任何樣本中都沒有足夠多的序列片段的基因應該從下游分析中過濾掉。這樣做的原因有好幾個逢防。 從生物學的角度來看叶沛，在任何條件下的表達水平都不具有生物學意義的基因都不值得關注，因此最好忽略忘朝。 從統(tǒng)計學的角度來看灰署，去除低表達計數(shù)基因使數(shù)據(jù)中的均值 - 方差關系可以得到更精確的估計，并且還減少了在觀察差異表達的下游分析中需要進行的統(tǒng)計檢驗的數(shù)量局嘁「然”
還有這篇：
單細胞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。現(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是歸一化的一種方法或南，代表每百萬堿基中每個轉錄本的count值。
注意：這個算法只是校正文庫差異艾君，而沒有校正基因長度差異采够。要注意我們分析的目的就是：比較一個基因在不同細胞的表達量差異，而不是考慮一個樣本中不同兩個基因的差異冰垄，因為"沒有兩片相同的樹葉”這個差異是正常的蹬癌。但是同一個基因由于某種條件發(fā)生了改變，背后的生物學意義是更值得探索的虹茶。
用起來很簡單逝薪，有現(xiàn)成的函數(shù)cpm() ，然后我們再用log將數(shù)據(jù)降個維度蝴罪，但保持原有數(shù)據(jù)形狀不變：log2(edgeR::cpm(dat)+1) 董济。
意思就是：cpm需要除以測序總reads數(shù)，而這個值作為分母會導致結果千差萬別要门，有的特別大有的很小虏肾。為了后面可視化不受極值的影響，用log轉換一下可以將數(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
# 和標準公式比較看看炒瘟，結果完全相同
> 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")
