日常瞎掰
??目前侧啼，普通RNA-seq測(cè)序已經(jīng)成為科研中一件非常稀松平常的事牛柒，這也是得益于其物美價(jià)廉的基礎(chǔ)。雖然已經(jīng)有不少成熟的RNA-seq分析流程存在痊乾，但并沒(méi)有形成一個(gè)統(tǒng)一的標(biāo)準(zhǔn)答案皮壁。畢竟，條條大路通羅馬哪审，走哪一條都可以得到想要的結(jié)果蛾魄，這也無(wú)傷大雅。就拿我們今天要說(shuō)的RNA-seq的標(biāo)準(zhǔn)化方法來(lái)說(shuō)就不止一種湿滓，可供選擇的方式就有CPM滴须、TPM、FPKM三種叽奥。當(dāng)然扔水，今天我們并不是要討論哪一種方式更好，而是想跟大家分享一個(gè)好用的工具來(lái)快速進(jìn)行標(biāo)準(zhǔn)化朝氓。rnanorm是一個(gè)基于python開(kāi)發(fā)的命令行工具魔市，可以很方便地來(lái)做以上三種方式的標(biāo)準(zhǔn)化操作。
安裝
??python包的安裝方式赵哲，很簡(jiǎn)單嘹狞，一行命令即可實(shí)現(xiàn)：
pip install rnanorm
使用
- 正常情況
??正常情況下，表達(dá)矩陣?yán)锩娴幕蛎ǔＪ褂肎TF里面的symobol誓竿、ensembl原始ID命名磅网，這樣最方便的方式就是我們給軟件兩個(gè)參數(shù)：原始矩陣和GTF文件即可：
head readcount.txt
FEATURE_ID sampleA1 sampleA2 sampleA3 sampleB1 sampleB2 sampleB3
ENSMUSG00000000001.1 4248 6254 3409 7387 8099 3804
ENSMUSG00000000003.15 0 0 5 0 0 0
ENSMUSG00000000028.15 2239 3068 1801 2581 3081 1466
ENSMUSG00000000031.16 0 0 0 3 0 0
ENSMUSG00000000037.17 57 67 79 1 2 0
ENSMUSG00000000049.11 0 0 0 0 6 0
rnanorm readcount.txt --annotation gencode.vM25.annotation.gtf --gene-id-attr gene_id --tpm-output sample.tpm.txt
??示例指定的是TPM方式，如果想要CPM筷屡、FPKM的結(jié)果可以分別使用參數(shù)--cpm-output涧偷、--fpkm-output來(lái)指定簸喂。這三個(gè)參數(shù)并不沖突，可以同時(shí)指定燎潮，得到三種標(biāo)準(zhǔn)化結(jié)果喻鳄。
- 其他情況
??有些時(shí)候，我們拿到的矩陣?yán)锩婊虻拿Q(chēng)可能不是原始的ID确封，如下面所示除呵，雖然知道使用的是ensembl，但并不是GTF里面原始的ID爪喘。因?yàn)槿コ薎D點(diǎn)后面的后綴颜曾，所以即使提供了GTF文件，軟件也是沒(méi)法完成矩陣標(biāo)準(zhǔn)化秉剑。
??TPM泛豪、FPKM標(biāo)準(zhǔn)化的操作，除了矩陣以外侦鹏，其實(shí)另外還需要每個(gè)基因的長(zhǎng)度信息诡曙，提供GTF文件就是為了獲取這個(gè)信息。對(duì)于提供基因長(zhǎng)度略水，也可以使用軟件--gene-lengths來(lái)提供价卤。像下面這樣基因ID與GTF中不完全一致的情況，要么把矩陣?yán)锩婊騃D轉(zhuǎn)化成與GTF一致后再用上面的命令做標(biāo)準(zhǔn)化渊涝，要么提供一個(gè)基因長(zhǎng)度的兩列表格使用下面的命令來(lái)標(biāo)準(zhǔn)化荠雕。
head readcount.txt
FEATURE_ID sampleA1 sampleA2 sampleA3 sampleB1 sampleB2 sampleB3
ENSMUSG00000000001 4248 6254 3409 7387 8099 3804
ENSMUSG00000000003 0 0 5 0 0 0
ENSMUSG00000000028 2239 3068 1801 2581 3081 1466
ENSMUSG00000000031 0 0 0 3 0 0
ENSMUSG00000000037 57 67 79 1 2 0
ENSMUSG00000000049 0 0 0 0 6 0
awk -F '[\t"]' -v OFS="\t" 'BEGIN{print"FEATURE_ID","GENE_LENGTHS"}/^chr/&&$3=="gene"{len=$5-$4+1;print substr($10,1,18),len}' gencode.vM25.annotation.gtf >mm10.genelen.vM25.txt
head mm10.genelen.vM25.txt
FEATURE_ID GENE_LENGTHS
ENSMUSG00000102693 1070
ENSMUSG00000064842 110
ENSMUSG00000051951 465598
ENSMUSG00000102851 480
ENSMUSG00000103377 2819
ENSMUSG00000104017 2233
rnanorm readcount.txt --gene-lengths=mm10.genelen.vM25.txt --tpm-output=readcount.tpm.txt
??如果想要的標(biāo)準(zhǔn)化方式是CPM，可以忽略--annotation或者--gene-lengths參數(shù)驶赏，直接提供原始矩陣即可。
結(jié)束語(yǔ)
??從上面的示例可以看出使用rnanorm做標(biāo)準(zhǔn)化還是很方便的既鞠，不過(guò)有一些注意事項(xiàng)還是需要指出的：表達(dá)矩陣的基因那一列的列名必須為FEATURE_ID煤傍；提供GTF文件時(shí)，需要使用--gene-id-attr參數(shù)來(lái)指定基因ID的屬性名嘱蛋，如示例矩陣中基因名為ensembl蚯姆，則其在GTF中對(duì)應(yīng)的ID屬性為gene_id，如果是symbol的話(huà)就要換成gene_name洒敏；提供基因長(zhǎng)度文件時(shí)龄恋，兩個(gè)列名需為FEATURE_ID、GENE_LENGTHS凶伙。當(dāng)然郭毕，CPM、TPM函荣、FPKM標(biāo)準(zhǔn)化并不難显押，知道了方法自己完全可以寫(xiě)個(gè)函數(shù)來(lái)計(jì)算扳肛，但既然有了現(xiàn)成好用的工具，直接使用豈不是更省事省力乘碑！為軟件作者口頭點(diǎn)個(gè)贊挖息！
往期回顧
rGREAT | 基因組區(qū)間功能富集
利用UCSC預(yù)測(cè)啟動(dòng)子序列可能結(jié)合的轉(zhuǎn)錄因子
Vision | scRNA細(xì)胞相似性 + 差異signature
HiC | contacts vs distance
hdWGCNA | 單細(xì)胞數(shù)據(jù)共表達(dá)網(wǎng)絡(luò)分析
