泛癌分析基因表達矩陣處理

Note：數(shù)據(jù)從UCSC Xena下載，一共33種癌癥，這里對于如何下載不做介紹。

1. 文件下載和處理

下載好的數(shù)據(jù)共33個文件，放到同一個文件夾。因為是gz格式的文件，所以全部解壓，并刪除原壓縮文件。
建議從命令提示行解壓并刪除原壓縮文件，命令如下：

$ 先用cd 命令進入文件所在路徑
del *.gz   # 刪除命令
$ 然后批量重命名，將tsv格式改為txt格式
ren *.tsv *.txt # 重命名

處理后文件如下圖：

33個表達矩陣文件

2. 將33個文件中Ensemble ID轉(zhuǎn)為基因名

這里使用python處理。建議安裝Anaconda，安裝教程不做介紹。
需要下載人類基因組gtf文件，可以從Ensemble官網(wǎng)下載。下載很慢，這里付上鏈接，下載后自行解壓。
鏈接：https://pan.baidu.com/s/14QjLGi-IZ7CKK3XEUqGRuA
提取碼：3oak

代碼如下：

import pandas as pd
import os
import sys
import re

gtf = sys.argv[1]
sourceDir = sys.argv[2]
targetDir = sys.argv[3]

# 構建Ensemble ID 與geneSymbol的字典
gene_dict = {}
with open(gtf, 'r') as GTF:
    lines = GTF.readlines()
    for line in lines:
        pattern = re.compile(r'.+gene_id "(.+?)";.+gene_name "(.+?)";.+gene_biotype "(.+?)".+')
        results = pattern.match(line)
        if results:
            gene_dict[results.group(1)] = results.group(2)


for dirName, subFolders, fileNames in os.walk(sourceDir):
    for file in fileNames:
        filepath = os.path.join(dirName, file)
        # 文件名中提取腫瘤類型
        filepattern = re.compile('.+?-(.+?)[.]')
        results = filepattern.match(file)
        tumortype = results.group(1)
        # id轉(zhuǎn)為gene symbol
        targetpath = os.path.join(targetDir, tumortype+".txt")
        with open(targetpath, 'w') as WF:
            with open(filepath, "r") as RF:
                line_1 = RF.readline()
                WF.write(line_1)
                lines = RF.readlines()
                for line in lines:
                    li_lst = line.split("\t")
                    ensemble_ID = li_lst[0].split(r".")[0]
                    if ensemble_ID in gene_dict.keys():
                        li_lst[0] = gene_dict[ensemble_ID]
                        WF.write("\t".join(li_lst))

代碼復制到txt文件中保存為.py格式即可，這里命名為id2symbol。如果不會用pycharm的話就從Anaconda中運行代碼。
如下圖：

Anaconda Prompt

代碼需要三個參數(shù)：基因組文件名（包含完整路徑）、表達矩陣文件所在路徑、轉(zhuǎn)換后的文件輸出路徑
從Anaconda Prompt中進入代碼所在的路徑，運行命令，命令格式如下：

python 代碼名稱 基因組文件名（包含完整路徑） 表達矩陣文件所在路徑  轉(zhuǎn)換后的文件輸出路徑

如下圖（示例）：

代碼運行

文件較多，運行時間較長。運行后文件以腫瘤類型命名，如下圖：

id轉(zhuǎn)換后的文件

3. 從33個文件中提取單基因表達數(shù)據(jù)

為了方面后續(xù)在R中的可視化，提取的數(shù)據(jù)最后合并為長格式的數(shù)據(jù)
代碼如下：

import pandas as pd
import os
import sys
import re
import numpy as np

gene = sys.argv[1]
sourceDir = sys.argv[2]
targetDir = sys.argv[3]

def addTissueType(df):
    for index, row in df.iterrows():
        sample_id = int(index.split('-')[3][0:2])
        if sample_id >= 10:
            df.loc[index,'tissueType'] = 'normal'
        elif sample_id < 10:
            df.loc[index,'tissueType'] = 'tumor'

allDF = pd.DataFrame(columns = (gene, 'tissueType', 'tumorType'))
for dirName, subFolders, fileNames in os.walk(sourceDir):
    for file in fileNames:
        tumorType = file.split('.')[0]
        filepath = os.path.join(dirName, file)
        df = pd.read_table(filepath, index_col=0)
        df = df.loc[gene].to_frame()
        addTissueType(df)
        df['tumorType'] = tumorType
        allDF = pd.concat([allDF, df])

targetpath = os.path.join(targetDir, 'result.txt')
allDF.to_csv(targetpath, sep='\t')

代碼的運行方式同上。代碼運行的命令格式為：

python  代碼名稱  你想要提取的基因名  id轉(zhuǎn)換后的文件所在路徑  結(jié)果輸出路徑

比如想提取的基因是HIF1A

運行示例

輸出的文件為result.txt,打開如下：

結(jié)果文件