2016年12月25日日曜日

リファレンスゲノムのバージョン変換 - イネ編


イネのデータベースのリファレンスゲノムがIRGSP build 4あるいはIRGSP  build 5になっている場合には、今使っている リファレンスゲノムIRGSP1.0に変換する必要があります

Samtoolsを使うと簡単!FASTAファイルから部分配列を抜き出すを参考にして、リファレンスゲノムを切り出します

リファレンスにインデックスを付与する
samtools faidx ref_IRGSP4.fasta

FASTAのエントリー名と位置を指定して部分配列を切り出す
samtools faidx ref_IRGSP4.fasta chr1:12345-12445
>chr1:12345-12445
ATGA.... 

抜き出した部分配列を用いて、blast検索を行います
手法は以下の論文を参考にしました

Ohyanagi H, Ebata T, Huang X, Gong H, Fujita M, Mochizuki T, Toyoda A,
Fujiyama A, Kaminuma E, Nakamura Y, Feng Q, Wang ZX, Han B, Kurata N.
OryzaGenome: Genome Diversity Database of Wild Oryza Species. Plant Cell Physiol.
2016 Jan;57(1):e1. doi: 10.1093/pcp/pcv171. PubMed PMID: 26578696; PubMed Central
PMCID: PMC4722174.

(i) for each SNP, extract the 201 bp flanking sequence (the SNP nucleotide itself and the flanking 100 bp nucleotides on both sides) on the original genome (IRGSP-build4.0)
(ii) search the counterpart of the 201 bp on the latest genome (Os-Nipponbare-Reference-IRGSP-1.0) with BLASTN homology search, allowing no indel and one mismatch at most, with 100% coverage; (iii) if there is more than one homologous region, discard the SNP
(iv) if the SNP was converted onto a different chromosome, discard the SNP.


blastdb作成
makeblastdb -dbtype nucl -hash_index -in IRGSP-1.0_genome.fasta

blast
blastn -db IRGSP-1.0_genome.fasta -query test.fasta -out test_blastn.out -outfmt 6 -max_target_seqs 1

オプションは
-outfmt 6 タブ区切りで出力
-out filename



2016年12月24日土曜日

Mapchart 2.30の使い方-データ準備編

Mapchart 2.30のインストールについては、インストール編をご利用ください。

MapChartの書式はMapChart example 4のファイルをを参考にすると使いやすいです。

  • 連鎖地図の書式
group B
TG608    0
TG189    2
CT63B   21
TG1B    33
TG234   41
TG14    51
CT24   124
一行目のgroup Bは染色体や連鎖群の名前などを記入します

1列目マーカー
2列目:位置(遺伝距離)

  • マーカーの書式
group A
C1    0.0 C1
C2    2.0 C2
C3    22.1 C3
C4    29.1 C4
C5    30.9 C5
C6      47.9 C6
C7    50.3 C7
C8     61.0 C8
C9     74.2 C9
C10      94.2 C10
C11    105.2 C11

1列目マーカー
2列目:位置(遺伝距離)
3列目:書式設定

3列目の書式設定は、例えば
B:bold
I:italic

C2:color code 2
です。color codeは数字を変えると色が変更になります。


  • QTL intervalsの書式
qtls ;start of qtls section of linkage group A 
Ca-14 1.3 7.9 12.3 18.6
      ;first QTL: name Ca-14
      ;           start of outer (1.3) and inner (7.9) intervals,
      ;           end of inner (12.3) and outer (18.6) intervals,
      ;           no formatting, displayed with default color (C1),
      ;           default fill and line style (solid) and default font style

Ca-14はQTLの名前
1.3 7.9 12.3 18.6はQTLの位置を示します。
   
vg-1 73.0 73.0 85.0 85.0  C3 F4 
start of outer intervalinner intervalを同じにすれば線なし

     ;third QTL: only inner interval shown, since it overlaps outer interval (73.0 - 85.0)
     ;           color number 3 (C3), 四角の色の変更 デフォルトは黒

     ;           fill style F4 (diagonal hatch) 模様を変更可能
     ;the first, second and third QTLs don't overlap each other and are shown in the same
     ;"QTL band" (in the same vertical line)



2016年11月15日火曜日

QTLの寄与率


QTLの寄与率をr/qtlで計算しています。

以下のページを参考にしています
How can I get the the percentage of phenotype variance explained by a QTL?
https://groups.google.com/forum/#!topic/rqtl-disc/wyHxuaBr-gA

fitqtl() を使って計算します。

formula="y~Q1+Q2+Q3"

out.fq <- fitqtl(cross, qtl=qtl, formula=formula,pheno=trait.id)

summary(out.fq)

$result.full
     df       SS        MS      LOD     %var Pvalue(Chi2)    Pvalue(F)
Model   3 1515.210 505.06993 8.056636 17.08753  4.37793e-08 6.030568e-08
Error 194 7352.131  37.89758       NA       NA           NA           NA
Total 197 8867.341        NA       NA       NA           NA           NA

$result.drop
     df Type III SS      LOD     %var  F value Pvalue(Chi2)    Pvalue(F)
2@23.4  1    560.4122 3.158383 6.319958 14.78754 0.0001368593 0.0001631004
6@34.5  1    474.0002 2.686255 5.345460 12.50740 0.0004361334 0.0005070081
8@46.7  1    358.9511 2.049507 4.048013  9.47161 0.0021249709 0.0023887012

QTLの寄与率は
%var shows QTL contributions(%).

2@23.4  -> 6.319958
6@34.5  -> 5.345460
8@46.7 -> 4.048013

QTLー座の分散の、表現型分散に対する比を QTLの寄与率としています

%var = percent of the phenotype variance explained by the QTL. 

2@23.4  ->  6.319958 =  560.4122/8867.341


QTL三座の寄与率は17.08%です



2016年10月25日火曜日

Mapchart 2.30の使い方-インストール編-

QTLの情報を連鎖地図に書き込みたいときに便利なツールをご紹介します。

Mapchart 2.30
Software for the graphical presentation of linkage maps and QTLs
https://www.wur.nl/en/show/Mapchart-2.30.htm

Join Mapからの出力でなくても、自前の連鎖地図とQTL情報を書き込み、作図することが可能です。
残念ながら、Windowsのみで使うことが可能です。

まずはDownload
のページからダウンロード可能です。
ソフトをダウンロードするだけでなく、license fileが必要です。




登録後入手したMapChart.licというファイルを、MapChartがインストールされているディレクトリに移動すれば、使用可能です。

つぎのデータ準備編へ

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

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現在の内容です。
更新などがある可能性がございますので、実行時にはご注意ください。

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

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で個体を指定するのは面倒ですが、ファイルを作るよりは簡単かもしれません


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


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

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

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]

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