2016年7月2日土曜日

plink formatをvcfに変換 - convert plink files to VCF

PLINK 1.90 betaを準備

(20180208追記)

ダウンロードはここから
https://www.cog-genomics.org/plink2/

必ずPLINK 1.90を使ってください。
# Convert binary files to vcf format. 

--bfileはPLINKのファイルがbinary fileの時に使います。


plink --file your_ped_map_input --recode vcf
plink --bfile your_bed_bim_fam_input --recode vcf
投稿者 時刻: 0 件のコメント:

2016年4月22日金曜日

AWS S3とフォルダを同期する


AWSコマンドラインインターフェースの利用

AWS Command Line Interface　OS Xで設定
以下のページを参考にしました
http://docs.aws.amazon.com/cli/latest/userguide/installing.html
Python softwareのpipを利用してインストールしますが、Python 2 version 2.6.5+ あるいは Python 3 version 3.3+が必要

$ sudo pip install awscli --ignore-installed six

AWS CLI の設定
aws configure コマンドを使って、セットアップ
以下のページを参考にしました。
$ aws configure
AWS Access Key ID [None]: AKIAI**********
AWS Secret Access Key [None]: wJalrXUtn*********
Default region name [None]: ap-northeast-1
Default output format [None]: json

Default region name [None]は以下のページから確認します。
Asia Pacific (Tokyo)なのでap-northeast-1を入力
http://docs.aws.amazon.com/general/latest/gr/rande.html#s3_region
Default output format は jsonを指定

LinuxからAWS S3へのアップロードは以下のページを参考にしました。
http://www.checksite.jp/amazon-s3-linux-filemanage/
S3上のフォルダの中を確認
$ aws s3 ls s3://バケット名/フォルダ名/

アップロードは「sync」コマンドを使って、ローカルフォルダと同期する形で
$ aws s3 sync test_folder/ s3://バケット名/フォルダ名/

「cp」コマンドを使って、ファイルのコピーも可能
$ aws s3 cp test.txt  s3://バケット名/フォルダ名/
crontabコマンドで定期的にフォルダの中身を更新することも可能のようですが、まだ設定していません...

この内容は2016.4.22現在の内容です。
更新などがある可能性がございますので、実行時にはご注意ください。

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

2016年4月2日土曜日

シェルスクリプトで文字列をランダムに並び替え


1)配列を入力
hkane$ array=(A1 A11 B54 C65 A23 A45 A33 A65 A78)

2) 配列をまずは配列を表示
全要素を改行区切りで配列に設定します
for item in "${array[@]}"; do echo "$item"; done
以下のように表示されます。
A1
A11
B54
C65
A23
A45
A33
A65
A78

3) ランダムに並べ替え
sortコマンドに-Rオプションがなかったため、テキストファイルをランダムに並べ替えるワンライナーを参考にして並び替え
while read x; do echo -e "$RANDOM\t$x"; done | sort -k1,1n | cut -f 2-
以下のように表示されます。
C65
A23
A33
A11
A1
A65
A78
B54
A45

4)元のarrayに戻します。
array=(`for item in "${array[@]}"; do echo "$item"; done | while read x; do echo -e "$RANDOM\t$x"; done | sort -k1,1n | cut -f 2-`)

5)arrayを表示させて確認します。
echo "${array[@]}"
A33 A45 A11 C65 A1 B54 A23 A78 A65


ランダムに配列を並べ替えて、5つだけ取り出す場合には
array=(`for item in "${array[@]}"; do echo "$item"; done | while read x; do echo -e "$RANDOM\t$x"; done | sort -k1,1n | cut -f 2-|head -n 5`)
と、head コマンドで指定可能です。

以下のページを参考にさせて頂きました。
ありがとうございます。

揮発性のメモ
http://d.hatena.ne.jp/iww/20130202/sort

配列の要素でループする
http://shellscript.sunone.me/array.html#配列の要素でループする:8177e845bec4a15c29be1b787e514e77

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

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 件のコメント:
登録: コメント (Atom)