上期(#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文件直接运行上面脚本即可.
更多原创精彩内容敬请关注生信杂谈：
