データ解析入門2018

データ解析入門2018

Windows Subsystem for Linux (WSL)のインストール

以下の点を確認する

  1. Windowsを最新版にアップデート (2018年9月だとバージョン1803であること! 確認方法: キーボードのWindowsキーを押しながら、[R]キーを押し、[ファイル名を指定して実行]画面が表示されるので、[名前]欄に「winver」と入力し、[OK]ボタンを選択します。)
  2. WSLのインストール 次のURL等を参照 https://qiita.com/yoshikingt/items/ab86411e6031459db805
  3. 日本語化: WSLの日本語化を参照

Hello World!

変数

if文

演算子

FASTQ, FASTAのリード数カウント

for文

配列

配列のソート

FASTA -> Tab区切り形式への変換

正規表現

複数ファイルの処理

Blast結果を発現情報にアノテーション

正規表現を用いた置換

VCFファイル整形

エスケープ処理

注意が必要な文字を含むテキストの処理

二次元配列

アンプリコンデータ解析

関数

GOツリー構築

絶対参照・相対参照

練習問題1

下記のようにサンプルごとに各遺伝子の発現量の表があるとする。この表はRNA-seqのリード数をそのまま表示したRAWデータだとする。これを補正してみる。

sample1sample2sample3sample4
gene1178199666569
gene2779047117
gene3116516216904
gene4222720471345

1.まずはサンプル毎にリード数が違うとサンプル間の比較ができないので、サンプル毎に合計リード数を算出する。(sum関数を使用する)

2.次に各サンプルの遺伝子ごとに次のようなRPM (Reads Per Million mapped reads)という値に変換した表を作成せよ。

( (各サンプルの遺伝子ごとのリード数) / (各サンプルの合計リード数) ) * 1000 * 1000

テキストファイルをExcelで読み込む

フィルター

練習問題2

4日目に行ったblastの情報を発現量データにアノテーションした結果ファイルを下記のコマンドでダウンロードする。

wget http://www.suikou.fs.a.u-tokyo.ac.jp/yosh_data/2018train/take.txt

1.フィルターを使用して、「pearl」を含む遺伝子の発現量のみ表示せよ。

2.さらに、その中で遺伝子発現量の平均値が高い順に表示せよ。


条件付き書式

練習問題3

練習問題2で使用したファイルに条件付き書式を使用して、発現量の大小が一目でわかるようにせよ。


関数、マクロ、VBA

練習問題4

練習問題2で使用したファイルにマクロを使用して、1~50行までの間、1行ずつ条件付き書式を使用して、発現量の大小が一目でわかるようにせよ。

回答例:take-macro.xlsm


練習問題5

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

FASTAファイル相同性検索blast, last
ゲノム(ショートリード)bwa, bowtie2
ゲノム(ロングリード)minimap2, minialign
RNA-seqtophat, hisat

ツールのインストール

上記のツールは比較的メジャーなため、実行ファイルをダウンロードするだけで良いときもあるが、Perl, Python, AWKなどのスクリプトファイルになっていることも多く、それらのインタープリタをインストールする必要がある。また、自分でソースファイルから実行ファイルをビルドする(作成する)必要があることもある。今回ツールとは別に必要な環境は下記のコマンドでインストール出来る。

sudo apt update
sudo apt install build-essential unzip python zlib1g-dev

ツールごとにインストール方法は異なるので、ツールのマニュアルをよく読んでインストールすること。

マニュアルはツールのHPに纏められていたり、ダウンロードしたファイルの中にREADME, MANUALといった名前のファイルに書いてあったりする。もしくは、ツールのオプションとして、「-h」、「–help」などをつけて実行したら出てくる説明文を読んだほうが早い場合もある。マニュアルを探す順番は、本家WEBサイト→ダウンロードしたファイルの中→とりあえず実行して表示されるメッセージ→本家以外のサイトをググる という順番が無難。ただ慣れないうちは、日本語で使い方が書いてあるページをググってみるのも良いかも。本家サイト以外では情報が古いことが多く、現在のバージョンと違いがある可能性を考慮すること。

ダウンロードしたファイルは基本的には圧縮されているため、フォーマットに適した解凍を行う必要がある。

ファイル名解凍方法
xxx.tar.gztar zvxf xxx.tar.gz
xxx.gzgzip -d xxx.gz
xxx.tar.bz2tar jvxf xxx.tar.bz2
xxx.bz2bzip2 -d xxx.bz2
xxx.tartar vxf xxx.tar
xxx.zipunzip 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-intronRNA-seqでイントロンを挟まない配列
rnaseq-intronRNA-seqでイントロンを挟む配列
p53_partialショウジョウバエのp53遺伝子イントロンがあります
nanoporeナノポアのデータ(低クオリティ、ロングリード)

大まかには次のキーワードに注意して挙動の違いを見て下さい。

  • 同じスコアで複数箇所にヒットした場合
  • 何回か実行してみて同じ結果が得られるか
  • グローバルアライメントとローカルアライメント
  • イントロンを含むリード(大きなdeletionがあるように見える)

出力ファイル形式

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

2018mapping例

Rのインストール

  1. R本体のインストール

    https://cran.r-project.org/

ライブラリのインストール

Rの本体だけではあまり大したことは出来ないので、研究者の人たちが作成したライブラリをインストールして使用することになる。

ライブラリのインストールは各ライブラリのHPの通りに行えばよいが、基本的には下記の2通りで行うことになる。

  1. CRANから(例:veganをインストールする場合)

    install.packages("vegan")

  2. Bioconductorから(例:DESeq2をインストールする場合)
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2", version = "3.8")

インストールが終わったら、いずれの方法でインストールしても、

library("DESeq2")

などとライブラリ名を打てばライブラリがロードされる。

Rの用途

データ解析におけるRの主な用途は、統計検定やグラフ作成になる。グラフ作成時には、どのようなグラフをRで描くことが可能か、下記のようなサイトで調べると良い。

R言語入門

R言語をざっくり読めるようになるには、下記のサイトなどが良さそう。

https://qiita.com/stkdev/items/6aba2c1db2fa056170ae

Rで変数の中身を理解するには

  • RStudioで変数名をクリック
  • class() 関数
  • str() 関数

①普通に作図:散布図の作成

サンプルデータ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つ選んで描いてみる

②ggplot2を使って作図:BOX AND SCATTER PLOT

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個だけ使用してヒートマップを描く

  • データ解析入門2018.txt
  • 最終更新: 2019/07/15 02:46
  • by 127.0.0.1