1.5日目練習問題で使用したexample.vcfファイルの10列目、11列目にはそれぞれ「sample1」、「sample2」の検体の変異情報が記載されている。
「:」区切りの最初の1項目目がジェノタイプ情報であり、下記の意味になる。 0/0 : 変異なし 0/1 : 片方(父親由来もしくは母親由来のどちらか)の染色体に変異あり(ヘテロ) 1/1 : 両方の染色体に変異あり(ホモ) ./. : デプスが不十分で判定不可能
ここで数字の1が2、3、4・・・などになることがあるが、それはVCFの5列目(ALT)の変異した配列の何番目になるかを示している。例:ALT→A,Cのとき、0/1であれば、Aの変異を1つ、0/2であればCの変異を1つ持つ。
sample1でデプスが不十分で判定不可能とされる変異はいくつか?sample2では?
#sample1
grep -v "^#" example.vcf|awk -F'\t' '{if($10~"^[.]/[.]"){cnt++}} END{print cnt}'
#sample2
grep -v "^#" example.vcf|awk -F'\t' '{if($11~"^[.]/[.]"){cnt++}} END{print cnt}'
また、DP(各検体のリード数)が4以下の変異数はsample1, 2それぞれ幾つか?
#sample1
grep -v "^#" example.vcf|awk -F'\t' '
{
n++;
split($10,arr,":");
if(arr[3]<=4){cnt++}
}
END{print n, cnt}
'
2.4日目で使用した、take.blastn.txtにはRで読み込むと意図しない結果になる文字(#)やほかのプログラムでエラーになりそうな文字が含まれている。次のそれぞれの文字が含まれている行数を出し、すべて「_」に置換せよ。
# ' ( )
awk -F'\t' '
{
if($0~"[#'"'"'()]"){cnt++};
gsub("[#'"'"'()]","_",$0);
print $0
}
END{print cnt}
' take.blastn.txt
3.とあるプログラムを使うためにtake.blastn.txtの2列目の遺伝子名を「“」で囲む必要が出てきた。2列目を「”」で囲んで出力せよ。
awk -F'\t' '
{
OFS="\t"; #$2を変更するともともとタブ区切りだったのがスペースに変換されてしまうので、タブのまま出力するのに必要
$2="\""$2"\"";
print $0
}
' take.blastn.txt