上期(#TCGA系列#TCGA基因/miRNA表達(dá)譜及臨床數(shù)據(jù)下載?)介紹了使用TCGA 的API下載腫瘤表達(dá)譜及臨床數(shù)據(jù),本期來處理上期下載的表達(dá)譜文件.還是以肝癌的miRNA表達(dá)譜為例.

我們上次已經(jīng)下載了373個cases的425個表達(dá)譜文件,每個樣本(case)的表達(dá)譜文件格式如下.

單個樣本miRNA表達(dá)譜

其他所有樣本的格式與上圖相同.每列依次是miRNA名稱,原始reads數(shù)目,歸一化reads數(shù)RPM,最后一列cross-mapped miRNA.

目錄結(jié)構(gòu)如下,都是file_ID/file_name的:

425個表達(dá)譜文件結(jié)構(gòu)

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

包含file_ID和file name的manifest文件

我們的目的是將425個表達(dá)譜文件合并成一個表達(dá)譜矩陣,并且以file_ID為列名,如結(jié)果是類似下面的:

表達(dá)量矩陣

shell腳本如下:

# 合并425個樣本的miRNA名及對應(yīng)表達(dá)量RPM,最終結(jié)果應(yīng)該是1882行miRNA和425列樣本表達(dá)量的矩陣文件,代碼如下:

# file_ID和file_name數(shù)組分別存儲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`)

# 數(shù)組file_path存儲文件路徑:
for((i=0;i<${#file_ID[@]};i++)){
    file_path[$i]="./"${file_ID[$i]}"/"${file_name[$i]}
    echo ${file_path[$i]}
}

# 使用awk二維數(shù)組進(jìn)行合并:
awk -v file_num=${#file_path[@]} '
    BEGIN{
        OFS="\t";
    }
    {
        # 每一個文件第一行是列名,而我們不需要合并列名,所以要NR>1
        # 然后以miRNA($1),文件ID(ARGIND),構(gòu)建值為表達(dá)量($2)二位數(shù)組a[miRNA][exp].
        if(FNR>1){a[$1][ARGIND]=$3;}
    }
    # 構(gòu)建了425個數(shù)組后進(jìn)行合并:
    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的表達(dá)量
                printf "%s\t",a[i][j];
                j=j+1;
            }
            print ""    #每一行加個換行
        }
    }' ${file_path[@]} >../miRNA_exp_matrix.txt

# 將file_ID添加到表達(dá)量矩陣中:
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添加到表達(dá)量矩陣中也可以使用以下代碼:
aaa=`echo miRNA ${file_ID[@]}|sed 's/ /\t/g' `
sed -i "1i $aaa" ../miRNA_exp_matrix.txt 

這個腳本運(yùn)算速度很快,2s左右.多樣本基因表達(dá)譜整合也是如此,只需下載所有的單個表達(dá)譜文件后替換manifest文件直接運(yùn)行上面腳本即可.

更多原創(chuàng)精彩內(nèi)容敬請關(guān)注生信雜談：