其他
前面几期已经完成了从TCGA表达谱/临床数据下载处理整合,本次对样本中的有annotations.txt的进行分析处理
在之前下载的表达文件中,可能已经有人注意到有些样本文件夹中还有个annotations.txt文件.
先看看这样的annotations文件有多少:
# 查看annotations.txt的文件的数量
ls ./*/annotations.txt|wc -l # 结果是23发现有23个.随便打开一个是这样的:
再打开一个是这样的:
重点注意notes部分,是对这个样本的说明.不过这样一个个点开看不方便,我们编个脚本合并下:
#将 23个annotations.txt合并到一个文件:
anno_files=(`ls ./*/annotations.txt`)
awk 'BEGIN{OFS="\t"}{if(ARGIND==1){print $0};if(ARGIND>1 && FNR>1){print $0} }' ${anno_files[@]} >../annotations_all.tsv这样就得到了一个有全部23个annotations.txt文件的合并文件annotations_all.tsv,前面几列分别是注释ID,提交者ID,cases_ID,所属范围,分类,创建时间,状态,注释信息.前几列都不重要,重要的是最后一列注释信息也就是notes,23个样本的注释信息如下:
这些notes大概包括以下几方面的注释:
属于肝胆混合样本,适合在胆管癌中研究,不适合在肝癌中研究.
在TCGA统计前已经经历过治疗.
先前有结肠癌 或肾细胞癌 或膀胱癌等癌症,接受或没有接受治疗,或者治疗史不明
其他…
通过上面的说明可以看出这23个样本都不是正常的肝癌样本,为了防止这23个特例在分析过程中产生不可控的因素,而且相对于整体样本数量425来说,去除这23个样本也不构成影响.所以我们去除表达谱中这23个样本.然后再合并为表达量矩阵.
之前的讲过的内容(TCGA基因/miRNA表达谱数据整合)已经可以生成miRNA表达谱矩阵文件.通过对之前的那个脚本修改,得到的下面这个脚本会跳过有annotaions的样本然后合并为表达量矩阵.
将这个脚本存为
merge_except_anno.sh.这个脚本需要三个参数, 第一个是manifest文件路径,第二个是'reads'或'rpm',第三个是输出文件路劲.
#! /bin/bash
#需要三个参数,$1是manifest文件路径,$2是'reads'或'rpm' ,$3是输出文件路劲.
manifest_path=$1
exp_patten=$2
out_path=$3
file_ID=(`awk '{if(NR>1)print $1}' $manifest_path`)
file_name=(`awk '{if(NR>1)print $2}' $manifest_path`)
anno_files_ID=(`ls ./*/annotations.txt|awk 'BEGIN{FS="/"}{print $2}' -`)
# 排除anno_files_ID后,使用数组file_path存储文件路径:
for((i=0;i<${#file_ID[@]};i++)){
if [[ "${anno_files_ID[@]}" =~ "${file_ID[$i]}" ]]
then
echo "跳过注释文件:"${file_ID[$i]}
continue
fi
file_path[$i]="./"${file_ID[$i]}"/"${file_name[$i]}
}
echo "去除注释文件还有"${#file_path[@]}"个文件"
# 使用awk二维数组进行合并:
awk -v file_num=${#file_path[@]} -v exp_patten=$exp_patten '
BEGIN{
OFS="\t";
}
{
# 每一个文件第一行是列名,而我们不需要合并列名,所以要NR>1
# 然后以miRNA($1),文件ID(ARGIND),构建值为表达量($2)二位数组a[miRNA][exp].
if(FNR>1 && exp_patten=="rpm"){a[$1][ARGIND]=$3;}
if(FNR>1 && exp_patten=="reads"){a[$1][ARGIND]=$2;}
}
# 构建了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[@]} > $out_path
#添加首行列名:
aaa=`echo ${file_path[@]}|sed 's/ /\n/g'|awk 'BEGIN{FS="/";ORS="\t"}{print $2}'|awk 'BEGIN{printf "%s\t","miRNA"}{print}' - `
sed -i "1i $aaa" $out_path
echo "成功生成去除注释的$exp_patten模式的表达量矩阵文件:"$out_path生成
reads count作为表达量值的表达量矩阵文件exp_matrix_reads.tsv 的脚本:
../merge_except_anno.sh ../gdc_manifest.2017-05-26T16-02-11.963011.tsv "reads" ../exp_matrix_reads.tsv生成
reads_per_million作为表达量值的表达量矩阵文件exp_matrix_rpm.tsv 的脚本:
../merge_except_anno.sh ../gdc_manifest.2017-05-26T16-02-11.963011.tsv "rpm" ../exp_matrix_rpm.tsv本期查看并去除了有annotations的样本,并将剩余样本进行了整合,得到了最终的两个表达量矩阵文件,这为之后的差异表达分析及生存分析做好了准备.
更多原创精彩内容敬请关注生信杂谈：