hkaneのLab Notebook

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 intervalとinner 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

http://www.wur.nl/en/Expertise-Services/Research-Institutes/plant-research/Download-MapChart.htm

のページからダウンロード可能です。

ソフトをダウンロードするだけでなく、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

以下のように表示されます。

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

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