練習問題7
wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/sample1.fasta.gz wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/sample2.fasta.gz wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/sample3.fasta.gz wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/sample4.fasta.gz gzip -d sample*.fasta.gz
として練習用のVCFファイルをダウンロードして解凍する。 上記ファイルは抗体可変領域のアンプリコンシーケンスデータである。
3．サンプル1～4のファイルについて配列の出現回数を求め、下記のように統合された出現頻度一覧表を作成せよ。
回答例：
awk -F'\t' '
BEGIN{
ORS="";
}
{
if(FILENAME==ARGV[1]){
if(FNR%2==0){
data[1][$0]=data[1][$0]+1;
seq[$0]=1;
}
}
if(FILENAME==ARGV[2]){
if(FNR%2==0){
data[2][$0]=data[2][$0]+1;
seq[$0]=1;
}
}
if(FILENAME==ARGV[3]){
if(FNR%2==0){
data[3][$0]=data[3][$0]+1;
seq[$0]=1;
}
}
if(FILENAME==ARGV[4]){
if(FNR%2==0){
data[4][$0]=data[4][$0]+1;
seq[$0]=1;
}
}
}
END{
print "id\tsample1\tsample2\tsample3\tsample4\n";
for(i in seq){
print i;
for(j=1;j<=4;j=j+1){
print "\t"data[j][i]+0;
}
print "\n";
}
}' sample1.fasta sample2.fasta sample3.fasta sample4.fasta
より簡潔な回答例：
awk -F'\t' '
BEGIN{
ORS="";
}
{
for(i=1;i<=length(ARGV)-1;i++){ #ARGV[0]はプログラム名(awk)が入っているため、-1しないと引数の数にはならない
if(FILENAME==ARGV[i]){
if(FNR%2==0){
data[i][$0]=data[i][$0]+1;
seq[$0]=1;
}
}
}
}
END{
#ヘッダー行の出力
print "id";
for(j=1;j<=length(ARGV)-1;j=j+1){
print "\t"ARGV[j];
}
print "\n";
#数値データの出力
for(i in seq){
print i;
for(j=1;j<=length(ARGV)-1;j=j+1){
print "\t"data[j][i]+0;
}
print "\n";
}
}' sample*.fasta
より簡潔な回答例2：
awk -F'\t' '
BEGIN{
ORS="";
}
{
if(FNR%2==0){
data[FILENAME][$0]=data[FILENAME][$0]+1;
seq[$0]=1;
}
}
END{
#ヘッダー行の出力
print "id";
for(j=1;j<=length(ARGV)-1;j=j+1){
print "\t"ARGV[j];
}
print "\n";
#数値データの出力
for(i in seq){
print i;
for(j=1;j<=length(ARGV)-1;j=j+1){
print "\t"data[ARGV[j]][i]+0;
}
print "\n";
}
}' sample*.fasta
4．連結した一覧表の各サンプル当たりのリード数が1万になるように標準化せよ
5．補正されたリード数の平均値が多い順にテーブルをソートせよ。
回答例
awk -F'\t' '
BEGIN{
ORS="";
}
{
if(FNR%2==0){
data[FILENAME][$0]=data[FILENAME][$0]+1;
seq[$0]=1;
}
}
END{
#まずサンプルごとのリード数を算出する
for(j=1;j<=length(ARGV)-1;j=j+1){
for(i in data[ARGV[j]]){
total[j]=total[j]+data[ARGV[j]][i];
}
}
#サンプルごとに合計1万リードになるように補正
for(j=1;j<=length(ARGV)-1;j=j+1){
for(i in data[ARGV[j]]){
data[ARGV[j]][i]=data[ARGV[j]][i]*10000/total[j];
#ついでに全サンプル合計した配列ごとのカウント値を計算しておく
cnt[i]=cnt[i]+data[ARGV[j]][i];
}
}
#ヘッダー行の出力
print "id";
for(j=1;j<=length(ARGV)-1;j=j+1){
print "\t"ARGV[j];
}
print "\n";
#数値データの出力
PROCINFO["sorted_in"]="@val_num_desc";
for(i in cnt){
print i;
for(j=1;j<=length(ARGV)-1;j=j+1){
print "\t"data[ARGV[j]][i]+0;
}
print "\n";
}
}' sample*.fasta