#TCGA系列#TCGA基因/miRNA表达谱数据整合

上期(#TCGA系列#TCGA基因/miRNA表达谱及临床数据下载​)介绍了使用TCGA 的API下载肿瘤表达谱及临床数据,本期来处理上期下载的表达谱文件.还是以肝癌的miRNA表达谱为例.

我们上次已经下载了373个cases的425个表达谱文件,每个样本(case)的表达谱文件格式如下.

单个样本miRNA表达谱

其他所有样本的格式与上图相同.每列依次是miRNA名称,原始reads数目,归一化reads数RPM,最后一列cross-mapped miRNA.

目录结构如下,都是file_ID/file_name的:

425个表达谱文件结构

file_ID和file_name在上期下载的manifest中有,manifest文件如下:

包含file_ID和file name的manifest文件

我们的目的是将425个表达谱文件合并成一个表达谱矩阵,并且以file_ID为列名,如结果是类似下面的:

表达量矩阵

shell脚本如下:

# 合并425个样本的miRNA名及对应表达量RPM,最终结果应该是1882行miRNA和425列样本表达量的矩阵文件,代码如下:

# file_ID和file_name数组分别存储file ID和file name
bash
file_ID=(`awk '{if(NR>1)print $1}' ../gdc_manifest.2017-05-26T16-02-11.963011.tsv`)
file_name=(`awk '{if(NR>1)print $2}' ../gdc_manifest.2017-05-26T16-02-11.963011.tsv`)

# 数组file_path存储文件路径:
for((i=0;i<${#file_ID[@]};i++)){
    file_path[$i]="./"${file_ID[$i]}"/"${file_name[$i]}
    echo ${file_path[$i]}
}

# 使用awk二维数组进行合并:
awk -v file_num=${#file_path[@]} '
    BEGIN{
        OFS="\t";
    }
    {
        # 每一个文件第一行是列名,而我们不需要合并列名,所以要NR>1
        # 然后以miRNA($1),文件ID(ARGIND),构建值为表达量($2)二位数组a[miRNA][exp].
        if(FNR>1){a[$1][ARGIND]=$3;}
    }
    # 构建了425个数组后进行合并:
    END{
        for(i in a){    # 一维是miRNA,所以i就是miRNA
            printf "%s\t",i     #输出miRNA
            j=1;        # 为了不改变文件顺序所以使用渐加的方式循环
            while(j<file_num+1){        #循环输出每个样本中miRNA的表达量
                printf "%s\t",a[i][j];
                j=j+1;
            }
            print ""    #每一行加个换行
        }
    }' ${file_path[@]} >../miRNA_exp_matrix.txt

# 将file_ID添加到表达量矩阵中:
echo miRNA ${file_ID[@]}|sed 's/ /\t/g'|awk '{if(NR==FNR)print;if(NR>FNR)print}' -  ../miRNA_exp_matrix.txt >../miRNA_exp_matrix_tmp.txt
cp ../miRNA_exp_matrix_tmp.txt ../miRNA_exp_matrix.txt
#删除临时文件:
rm ../miRNA_exp_matrix_tmp.txt

# 将file_ID添加到表达量矩阵中也可以使用以下代码:
aaa=`echo miRNA ${file_ID[@]}|sed 's/ /\t/g' `
sed -i "1i $aaa" ../miRNA_exp_matrix.txt 

这个脚本运算速度很快,2s左右.多样本基因表达谱整合也是如此,只需下载所有的单个表达谱文件后替换manifest文件直接运行上面脚本即可.

更多原创精彩内容敬请关注生信杂谈：