2016年3月24日木曜日

vcf fileから複数個体を取り出す


vcftoolsを使うと、指定した個体のvcf fileを作成することが出来ます

 まずは個体を指定するためのファイルを作成
取り出す個体をtxt fileに書き込みます

 たとえば
A1
A2
A3

のように、改行しながら個体名を記入していきます。

vcftools --gzvcf test.vcf.gz --keep <filename> 

 で、系統名を記入したファイルを指定します。

数個体だけ取り出す場合は、ファイルを作成する必要は有りません。

 vcftools --gzvcf test.vcf --indv A1 --indv A2 --indv A3

 何度も--indvで個体を指定するのは面倒ですが、ファイルを作るよりは簡単かもしれません
投稿者 時刻: 0 件のコメント:
ラベル: ,

2016年2月8日月曜日

vcf fileをGT(genotype)の情報のみを取り出す

vcf fileからGT(genotype)の情報のみを取り出し、

 hom_RR -> 0
het_RA -> 1
hom_AA -> 2

 の値に変換するためのスクリプト
あらかじめ、biallelicにフィルタリングしてあります。
vcftools_0.1.13


### GT のみを抽出

vcftools　--gzvcf　test.vcf.gz --extract-FORMAT-info GT --out test

test.GT.FORMATというファイルが作成されます。
このファイルにはGT(genotype)情報のみが記載されています。

scoreに変換して、tsv fileで出力

sed "s/1|1/2/g" test.GT.FORMAT | sed "s/0|1/1/g" | sed "s/1|0/1/g" | sed "s/0|0/0/g" > test_score.tsv


sed "s/1[/]1/B/g" test.GT.FORMAT | sed "s/0[/]1/H/g" | sed "s/1[/]0/H/g" | sed "s/0[/]0/A/g" | sed "s/[.][/][.]/NA/g" > test_score.tsv


間違いなどがありましたら、お知らせいただければ助かります。

投稿者 時刻: 0 件のコメント:
ラベル: , ,

2016年1月25日月曜日

vcftoolsでフィルタリング

様々な欠測率でフィルタリングする時のスクリプト
あらかじめリストを作成しておき、そのリストの値を代入する時のシェルスクリプト
vcftools_0.1.13
備忘録として

DIR=/Volumes/vcf
file=all_chr_snp
vcftools=/Volumes/vcftools_0.1.13/bin/vcftools

cd $DIR
mkdir filtered_vcf
missing=("0.2" "0.4" "0.6" "0.8")
for ((i=0; i< ${#missing[@]}; ++i))
do
$vcftools --gzvcf ${DIR}/${file}.vcf.gz --min-alleles 2 --max-alleles 2 --minDP 3 --maxDP 100 --remove-indels --max-missing ${missing[$i]} --minQ 20 --recode --out ${DIR}/filtered_vcf/${file}_MS${missing[$i]}_bi_MQ20_DP3-100
done

投稿者 時刻: 0 件のコメント:

2016年1月14日木曜日

連鎖地図のマーカーの物理位置を調べる

論文に書かれている連鎖地図上のQTLの位置を知りたいということがあります。その物理位置を推定するには、まずは連鎖地図のマーカーの物理位置を知る必要があります。

イネやソルガムのQTLはGrameneで調べることが多いです。

Gramene QTL Database

たとえば、ソルガムの第8染色体上の出穂関連QTL AQGN050を調べる場合、Associated Markersの情報から、SSRの情報を得ます。


と書かれています。

この配列を超絶高速ゲノム配列検索し、物理位置を調べます。primerのような短い配列は、blastでヒットしないことも多いため、こちらで検索しています。
大量のプライマーの位置情報を得るときには、GoogleSpreadSheet を利用します。
リンク先のシートを参考にして頂ければ、位置情報だけでなく、PCRのサイズも調べることが可能です。




これで大まかな物理位置を調べることが出来ました。

SSRの場合はこのように比較的簡単なのですが、RFLPマーカーの場合は難しくなります。
maizeのRFLPマーカーを利用している場合、MaizeGDBのデータベースから配列を調べることが可能です。
例えば、umc107というRFLPマーカーの場合、このマーカーをキーワードに検索を行うと



locusの情報が出てきます。

Probes:

の部分には p-umc107 (via RFLP Hybridization)へのリンクがあります。


 このページのリンク先を見るとSTS化された配列を入手できました。

この配列データから、このRFLPマーカーの物理位置を調べることが可能です。

他にもRFLPマーカーの配列情報を入手する方法としてはNCBIのデータベースでマーカー名を入力して検索を行うことが考えられます。

先ほどのumc107をキーワードに検索すると、Geneでヒットします。二つの染色体がヒットしますが、今回は第1染色体の情報を得たいのでumc107aの配列を調べることにします。




Name/Gene IDDescriptionLocationAliases
ID: 100282170		uncharacterized LOC100282170 [Zea mays]Chromosome 5, NC_024463.1 (10941609..10947256, complement)ZEAMMB73_755942, GRMZM2G090172
ID: 100192908

uncharacterized LOC100192908 [Zea mays]

連鎖地図のマーカーの物理位置を調べることは、このように手作業で検索することが多いです。
もう少し効率の良い方法をご存知でしたら、教えて頂きたいと思います。



投稿者 時刻: 0 件のコメント:

2015年10月17日土曜日

Tasselでvcf fileが読み込めない

Tassel　5.0　 を使っているのですが、どうしても読み込めないvcf fileがありました。

http://www.maizegenetics.net/#!tassel/c17q9

マニュアルをみると、vcf formatとして、tasselでは

http://www .1000genomes.org/wiki/analysis/variant-call-format/vcf-variant-call-format-version-42
https://cseweb.ucsd.edu/classes/fa14/cse182-a/notes/VCFv4.2.pdf
を採用しているようです。

読み込むことの出来るvcf fileでは
##fileformat=VCFv4.2

読み込むことが出来ないvcf fileでは
##fileformat=VCFv4.1
##fileformat=VCFv4.0
となっており、これが原因で読み込むことが出来ないかもしれません。

【追記】
本来ならばVCF fileをそのまま読み込めるはずですが、うまく行かない場合はvcftoolsを用いてvcf fileをplink形式などに変換します。

http://vcftools.sourceforge.net/man_latest.html#COMPARISON%20OPTIONS

OUTPUT OTHER FORMATS
--plink のオプションを使います。

 変換したplink fileを用いれば、無事に解析が出来ました。


投稿者 時刻: 0 件のコメント:

2015年10月14日水曜日

遺伝研スパコンからのデータ転送

スパコンへのデータ転送

  • ターミナルを使用して、localからスパコンへ
/Volumes/xxxx/xxx.shというファイルを転送する場合、

scp /Volumes/xxxx/xxx.sh username@gw.ddbj.nig.ac.jp:/home/username/xxxxx/

パスワードを要求されるので、入力するとデータが指定したパスに転送される


  • ターミナルを使用して、スパコンからlocalへ
sftp username@gw.ddbj.nig.ac.jp
を入力すると、
Connected to gw.ddbj.nig.ac.jp.
が表示される。
sftp> get /home/username/xxxxx/xxxxx.sam


投稿者 時刻: 0 件のコメント:

2015年10月6日火曜日

ubuntuのjavaをアップデート

AWSを使っていて、ubuntuのjavaのバージョンが古かったのでアップデート

以下のページを参考にさせて頂きました。

http://www.yottanote.com/contents/linux/ubuntu/oraclejava_update.html


ubuntu@********:~$ java -version
java version "1.7.0_79"

もともとインストールされているjavaはjava version "1.7.0_79"


sudo add-apt-repository ppa:webupd8team/java
sudo apt-get update
sudo apt-get install oracle-java8-installer

再度バージョンを確認

ubuntu@********:~$ java -version
java version "1.8.0_60"


投稿者 時刻: 0 件のコメント:
登録: 投稿 (Atom)