泛癌分析基因表达矩阵处理

Note：数据从UCSC Xena下载，一共33种癌症，这里对于如何下载不做介绍。

1. 文件下载和处理

下载好的数据共33个文件，放到同一个文件夹。因为是gz格式的文件，所以全部解压，并删除原压缩文件。
建议从命令提示行解压并删除原压缩文件，命令如下：

$ 先用cd 命令进入文件所在路径
del *.gz   # 删除命令
$ 然后批量重命名，将tsv格式改为txt格式
ren *.tsv *.txt # 重命名

处理后文件如下图：

33个表达矩阵文件

2. 将33个文件中Ensemble ID转为基因名

这里使用python处理。建议安装Anaconda，安装教程不做介绍。
需要下载人类基因组gtf文件，可以从Ensemble官网下载。下载很慢，这里付上链接，下载后自行解压。
链接：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转为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

代码需要三个参数：基因组文件名（包含完整路径）、表达矩阵文件所在路径、转换后的文件输出路径
从Anaconda Prompt中进入代码所在的路径，运行命令，命令格式如下：

python 代码名称 基因组文件名（包含完整路径） 表达矩阵文件所在路径  转换后的文件输出路径

如下图（示例）：

代码运行

文件较多，运行时间较长。运行后文件以肿瘤类型命名，如下图：

id转换后的文件

3. 从33个文件中提取单基因表达数据

为了方面后续在R中的可视化，提取的数据最后合并为长格式的数据
代码如下：

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转换后的文件所在路径  结果输出路径

比如想提取的基因是HIF1A

运行示例

输出的文件为result.txt,打开如下：

结果文件