前面幾期已經(jīng)完成了從TCGA表達譜/臨床數(shù)據(jù)下載處理整合,本次對樣本中的有annotations.txt的進行分析處理
在之前下載的表達文件中,可能已經(jīng)有人注意到有些樣本文件夾中還有個annotations.txt文件.
先看看這樣的annotations文件有多少:
# 查看annotations.txt的文件的數(shù)量
ls ./*/annotations.txt|wc -l # 結果是23
發(fā)現(xiàn)有23個.隨便打開一個是這樣的:
annotations.txt文件
再打開一個是這樣的:
![annotations.txt文件]](http://upload-images.jianshu.io/upload_images/903467-5c53537dd01057f6.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
重點注意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,所屬范圍,分類,創(chuàng)建時間,狀態(tài),注釋信息.前幾列都不重要,重要的是最后一列注釋信息也就是notes,23個樣本的注釋信息如下:
23個annotations.txt的注釋
這些notes大概包括以下幾方面的注釋:
- 屬于肝膽混合樣本,適合在膽管癌中研究,不適合在肝癌中研究.
- 在TCGA統(tǒng)計前已經(jīng)經(jīng)歷過治療.
- 先前有結腸癌 或腎細胞癌 或膀胱癌等癌癥,接受或沒有接受治療,或者治療史不明
- 其他...
通過上面的說明可以看出這23個樣本都不是正常的肝癌樣本,為了防止這23個特例在分析過程中產(chǎn)生不可控的因素,而且相對于整體樣本數(shù)量425來說,去除這23個樣本也不構成影響.所以我們去除表達譜中這23個樣本.然后再合并為表達量矩陣.
之前的講過的內容(TCGA基因/miRNA表達譜數(shù)據(jù)整合?)已經(jīng)可以生成miRNA表達譜矩陣文件.通過對之前的那個腳本修改,得到的下面這個腳本會跳過有annotaions的樣本然后合并為表達量矩陣.
將這個腳本存為
merge_except_anno.sh.這個腳本需要三個參數(shù), 第一個是manifest文件路徑,第二個是'reads'或'rpm',第三個是輸出文件路勁.
#! /bin/bash
#需要三個參數(shù),$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后,使用數(shù)組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二維數(shù)組進行合并:
awk -v file_num=${#file_path[@]} -v exp_patten=$exp_patten '
BEGIN{
OFS="\t";
}
{
# 每一個文件第一行是列名,而我們不需要合并列名,所以要NR>1
# 然后以miRNA($1),文件ID(ARGIND),構建值為表達量($2)二位數(shù)組a[miRNA][exp].
if(FNR>1 && exp_patten=="rpm"){a[$1][ARGIND]=$3;}
if(FNR>1 && exp_patten=="reads"){a[$1][ARGIND]=$2;}
}
# 構建了425個數(shù)組后進行合并:
END{
for(i in a){ # 一維是miRNA,所以i就是miRNA
printf "%s\t",i #輸出miRNA
j=1; # 為了不改變文件順序所以使用漸加的方式循環(huán)
while(j<file_num+1){ #循環(huán)輸出每個樣本中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的樣本,并將剩余樣本進行了整合,得到了最終的兩個表達量矩陣文件,這為之后的差異表達分析及生存分析做好了準備.
更多原創(chuàng)精彩內容敬請關注生信雜談：
