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

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

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

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

file_ID和file_name在上期下载的manifest中有,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

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

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