Windows Subsystem for Linux (WSL)のインストール
以下の点を確認する
下記のようにサンプルごとに各遺伝子の発現量の表があるとする。この表はRNA-seqのリード数をそのまま表示したRAWデータだとする。これを補正してみる。
| sample1 | sample2 | sample3 | sample4 | |
| gene1 | 178 | 199 | 666 | 569 |
| gene2 | 7 | 790 | 471 | 17 |
| gene3 | 116 | 516 | 216 | 904 |
| gene4 | 222 | 720 | 471 | 345 |
1.まずはサンプル毎にリード数が違うとサンプル間の比較ができないので、サンプル毎に合計リード数を算出する。(sum関数を使用する)
2.次に各サンプルの遺伝子ごとに次のようなRPM (Reads Per Million mapped reads)という値に変換した表を作成せよ。
( (各サンプルの遺伝子ごとのリード数) / (各サンプルの合計リード数) ) * 1000 * 1000
4日目に行ったblastの情報を発現量データにアノテーションした結果ファイルを下記のコマンドでダウンロードする。
wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/take.txt
1.フィルターを使用して、「pearl」を含む遺伝子の発現量のみ表示せよ。
2.さらに、その中で遺伝子発現量の平均値が高い順に表示せよ。
練習問題2で使用したファイルに条件付き書式を使用して、発現量の大小が一目でわかるようにせよ。
VCFデータの整形
wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/sample.txt
上記ファイルは、黒色系統の生物(b)と黒色系統の中から生まれた白色の劣性突然変異体(w)のジェノタイプ情報である。
1.まずは上記ファイルをExcelで読み込む。
2.ホモ(1/1)、ヘテロ(0/1)、変異なし(0/0)、デプスが低くてジェノタイプ不明(./.)を分かりやすく色付けせよ。
3.白色の遺伝様式は劣性の遺伝様式を示すことが知られている。そのため、白色になる原因変異の近傍は、ホモ接合となっていることが想定される。そこで白色系統の個体でホモ(1/1および0/0)、黒色系統の個体でヘテロ(0/1)となる割合をそれぞれ計算せよ。(ヒント：countif関数を使うと便利)
4.候補変異がありそうな領域はどこか？
回答例：sample.xlsx
上記のツールは比較的メジャーなため、実行ファイルをダウンロードするだけで良いときもあるが、Perl, Python, AWKなどのスクリプトファイルになっていることも多く、それらのインタープリタをインストールする必要がある。また、自分でソースファイルから実行ファイルをビルドする(作成する)必要があることもある。今回ツールとは別に必要な環境は下記のコマンドでインストール出来る。
sudo apt update sudo apt install build-essential unzip python zlib1g-dev
ツールごとにインストール方法は異なるので、ツールのマニュアルをよく読んでインストールすること。
マニュアルはツールのHPに纏められていたり、ダウンロードしたファイルの中にREADME, MANUALといった名前のファイルに書いてあったりする。もしくは、ツールのオプションとして、「-h」、「–help」などをつけて実行したら出てくる説明文を読んだほうが早い場合もある。マニュアルを探す順番は、本家WEBサイト→ダウンロードしたファイルの中→とりあえず実行して表示されるメッセージ→本家以外のサイトをググる という順番が無難。ただ慣れないうちは、日本語で使い方が書いてあるページをググってみるのも良いかも。本家サイト以外では情報が古いことが多く、現在のバージョンと違いがある可能性を考慮すること。
ダウンロードしたファイルは基本的には圧縮されているため、フォーマットに適した解凍を行う必要がある。
| ファイル名 | 解凍方法 |
| xxx.tar.gz | tar zvxf xxx.tar.gz |
| xxx.gz | gzip -d xxx.gz |
| xxx.tar.bz2 | tar jvxf xxx.tar.bz2 |
| xxx.bz2 | bzip2 -d xxx.bz2 |
| xxx.tar | tar vxf xxx.tar |
| xxx.zip | unzip xxx.zip |
macの人は、自分でビルドすると苦労することが多いので、https://bioconda.github.io/index.html#set-up-channels の手順でbiocondaをインストールしてから、https://bioconda.github.io/recipes.html で希望のパッケージを見つけてインストールするのが無難。
マッピングツールを動かし、特徴をつかむ練習。 下記の配列データがショウジョウバエのゲノムのどこに由来する配列なのか、正しく結論を出してください。
FASTQ: read.fq
FASTA: read.fa (blast, last用)
ショウジョウバエのゲノム配列を下記からダウンロードする。
http://hgdownload.cse.ucsc.edu/goldenPath/dm6/bigZips/dm6.fa.gz
準備した配列は次のようになっています。
| 配列名 | ヒント | 見てほしいところ |
| unique-seq1 | ゲノム中に1か所だけ | |
| low-complexity-seq1 | かなりATリッチ | blastでは・・・ |
| centromere-seq1 | セントロメアの一部配列（かなり多コピー) | |
| unique-seq1-copy, low-complexity-seq1-copy, centromere-seq1-copy | 上記配列と全く同じ配列 | 配列は同じでもヒットする箇所は・・・ |
| rnaseq-no-intron | RNA-seqでイントロンを挟まない配列 | |
| rnaseq-intron | RNA-seqでイントロンを挟む配列 | |
| p53_partial | ショウジョウバエのp53遺伝子 | イントロンがあります |
| nanopore | ナノポアのデータ(低クオリティ、ロングリード) | |
大まかには次のキーワードに注意して挙動の違いを見て下さい。
SAMフォーマット
http://crusade1096.web.fc2.com/sam.html
tophat2ではそれをSAMを圧縮したBAMファイルが出力される。それを見るにはsamtoolsというツールが必要だが、
sudo apt install samtools
とやればインストール可能で、使用するときは、
samtools view tophat_out/accepted_hits.bam
などとすればよい。またFASTAファイル中の特定の範囲を抜き出して見たい場合、例えばchr2Lの1-10番目の塩基を表示させるには、
samtools faidx dm6.fa chr2L:1-10
Rの本体だけではあまり大したことは出来ないので、研究者の人たちが作成したライブラリをインストールして使用することになる。
ライブラリのインストールは各ライブラリのHPの通りに行えばよいが、基本的には下記の2通りで行うことになる。
install.packages("vegan")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2", version = "3.8")
インストールが終わったら、いずれの方法でインストールしても、
library("DESeq2")
などとライブラリ名を打てばライブラリがロードされる。
データ解析におけるRの主な用途は、統計検定やグラフ作成になる。グラフ作成時には、どのようなグラフをRで描くことが可能か、下記のようなサイトで調べると良い。
R言語をざっくり読めるようになるには、下記のサイトなどが良さそう。
サンプルデータmariom.tpm.txtは、0h vs 24hの比較で有意に発現上昇or減少した遺伝子の0h,24h,48h,・・・のデータ。まずは普通に散布図を描いてみる。
data=read.table("mariom.tpm.txt",header=T,row.names=1,sep="\t")
plot(x=data[,1],y=data[,2])
色々なパラメータ例：http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html
練習
・data[,1], data[,2]をlog2()で変換して対数の値で表示
・RNA-seqで良く見るMAプロットを描いてみる
MAプロットはX軸にA, Y軸にMを描いたもの
M = log2(サンプル1) – log2(サンプル2)
A = { log2(サンプル1) + log2(サンプル2) } / 2
・0h、24hのサンプルをそれぞれ1つ選んで描いてみる
http://www.r-graph-gallery.com/89-box-and-scatter-plot-with-ggplot2/
# Ggplot2 library
library(ggplot2)
# Data
names=c(rep("A", 80) , rep("B", 50) , rep("C", 70))
value=c( sample(2:5, 80 , replace=T) , sample(4:10, 50 , replace=T),
sample(1:7, 70 , replace=T) )
data=data.frame(names,value)
#Graph
qplot( x=names , y=value , data=data , fill=names) + geom_boxplot() + geom_jitter()
・①のmariom.tpm.txtでサンプルを3つ、遺伝子を100個選んで、そのboxplotを描いてみる。
data=read.table("mariom.tpm.txt",header=T,row.names=1,sep="\t")
library(Heatplus)
data=as.matrix(data)
corrdist = function(x) as.dist(1-cor(t(x)))
hclust.avl = function(x) hclust(x, method="average")
plot(regHeatmap(data, legend=2, dendrogram = list(clustfun=hclust.avl, distfun=corrdist)))
・遺伝子の数が多すぎてヒートマップが見えにくいので、20個だけ使用してヒートマップを描く