GNU datamash を使って転置

バイオインフォで扱うデータってcolumnが多い場合が多々あります。

例えば、10X Genomicsの公開データ(bam)をsamtoolsで見てみると。

kimoton@DESKTOP-BL78EM7:~$ samtools view http://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_possorted_genome_bam.bam | head -1
[knet_seek] SEEK_END is not supported for HTTP. Offset is unchanged.
[bam_header_read] EOF marker is absent. The input is probably truncated.
ST-K00126:314:HFYL2BBXX:7:2103:14996:4725       272     1       10001   1       3S95M   *       0       0       CCGTAGCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC      <)-7-)7<JF7<FA7---<<A-FFFF-A<<AA----FFAJJJAAAJJAFJJJF<JJJA---FFA-AFFFFJ<FJF-FFFA7FFFJ7JFFAJ7FAFF<A      NH:i:3  HI:i:3  AS:i:91 nM:i:1  RE:A:I  BC:Z:TCGTCACG   QT:Z:AAFFFJJJ   CR:Z:TTAACTCGTAGAAGGA   CY:Z:AAFFFJJJJJJJJJJJ   CB:Z:TTAACTCGTAGAAGGA-1 UR:Z:GTCCGGCGAC UY:Z:JJJJJJJJJJ UB:Z:GTCCGGCGAC RG:Z:pbmc8k:MissingLibrary:1:HFYL2BBXX:7

横長でとても見にくい。 データの内容がわかりづらい。。

こんな時はdatamash transposeに渡してやりましょう。

kimoton@DESKTOP-BL78EM7:~$ samtools view http://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_possorted_genome_bam.bam | head -1 | datamash transpose
[knet_seek] SEEK_END is not supported for HTTP. Offset is unchanged.
[bam_header_read] EOF marker is absent. The input is probably truncated.
ST-K00126:314:HFYL2BBXX:7:2103:14996:4725
272
1
10001
1
3S95M
*
0
0
CCGTAGCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC
<)-7-)7<JF7<FA7---<<A-FFFF-A<<AA----FFAJJJAAAJJAFJJJF<JJJA---FFA-AFFFFJ<FJF-FFFA7FFFJ7JFFAJ7FAFF<A
NH:i:3
HI:i:3
AS:i:91
nM:i:1
RE:A:I
BC:Z:TCGTCACG
QT:Z:AAFFFJJJ
CR:Z:TTAACTCGTAGAAGGA
CY:Z:AAFFFJJJJJJJJJJJ
CB:Z:TTAACTCGTAGAAGGA-1
UR:Z:GTCCGGCGAC
UY:Z:JJJJJJJJJJ
UB:Z:GTCCGGCGAC
RG:Z:pbmc8k:MissingLibrary:1:HFYL2BBXX:7

転置してるだけですけど、とっても見やすくなりました。

Rのt()と同じですね。同じですけど、linuxコマンドとしてパイプで繋げられるのはとっても便利。

逆順にもできます。

kimoton@DESKTOP-BL78EM7:~$ samtools view http://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_possorted_genome_bam.bam | head -1 | datamash reverse | datamash transpose
[knet_seek] SEEK_END is not supported for HTTP. Offset is unchanged.
[bam_header_read] EOF marker is absent. The input is probably truncated.
RG:Z:pbmc8k:MissingLibrary:1:HFYL2BBXX:7
UB:Z:GTCCGGCGAC
UY:Z:JJJJJJJJJJ
UR:Z:GTCCGGCGAC
CB:Z:TTAACTCGTAGAAGGA-1
CY:Z:AAFFFJJJJJJJJJJJ
CR:Z:TTAACTCGTAGAAGGA
QT:Z:AAFFFJJJ
BC:Z:TCGTCACG
RE:A:I
nM:i:1
AS:i:91
HI:i:3
NH:i:3
<)-7-)7<JF7<FA7---<<A-FFFF-A<<AA----FFAJJJAAAJJAFJJJF<JJJA---FFA-AFFFFJ<FJF-FFFA7FFFJ7JFFAJ7FAFF<A
CCGTAGCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC
0
0
*
3S95M
1
10001
1
272
ST-K00126:314:HFYL2BBXX:7:2103:14996:4725

インストール法

Ubuntu

sudo apt-get install datamash

RHEL

wget http://files.housegordon.org/datamash/bin/datamash-1.0.6-1.el6.x86_64.rpm
sudo rpm -i datamash-1.0.6-1.el6.x86_64.rpm

MacOS

brew install datamash

参考

datamash - GNU Project - Free Software Foundation

R 3.5.0 へのアップデート(PPA利用)

きっかけ

最近のRパッケージ(今回はsinglecell解析に使うパッケージ)がR 3.5.0でないと動かない。。

f:id:kimoppy126:20180614091358p:plain

Bioconductor - SingleCellExperiment (development version)

R 3.5.0入れたdocker imageを動かそうとも考えたけど、データのマウントとか面倒そう。。 WSLのシステム中に入れちゃおう。

PPA(パーソナル・パッケージ・アーカイブ)を使用してインストールします。 PPAについてはこちらを参考にしてください。

要するにこういうことです。

PPAはUbuntuユーザーのチームや個人がそれぞれ管理している非公式のApp Storeのようなもので、Ubuntuの公式レポジトリからはダンロードできないソフトウェアや最新のバージョンのソフトウェアを手に入れることができます。
参考:UbuntuのPPAて何? [Linuxの使い方] All About

インストール

1. システム内のパッケージを削除

apt-getを用いてインストールした R(r-cran-* とつくもの) の削除を削除する。

対象となるファイルのリストアップ。

dpkg -l | grep r-cran-

apt-getで削除

sudo apt-get remove r-cran-*

2. Michael Rutter's PPAの追加

aptレポジトリにrutterさんのPPAを追加します。 忘れずにapt-get updateしてください。

sudo add-apt-repository ppa:marutter/rrutter3.5
sudo apt-get update

3. Rのアップグレード

取り込んだPPAからインストールします。

sudo apt install r-api-3.5

4. ライブラリの再インストール

R 3.4で使用していたライブラリを再インストールします。

installed <- rownames(installed.packages())
pkgs <- dir("~/R/x86_64-pc-linux-gnu-library/3.4")
length(installed)
[1] 14
length(pkgs)
[1] 280

R 3.4.4では280のライブラリがインストールされていました。 新しくインストールしたR3.5.0には初期状態で14のライブラリがインストールされています。

new <- setdiff(pkgs, installed)
length(new)
[1] 280

これら2つで共通していないものを抽出し、インストールします(今回は280個すべて)。 Ncpusで使用するCPU数を指定できます。

install.packages(new, Ncpus = 6)

結果

いくつかインストールに失敗してましたが、、

 警告メッセージ:
1:  パッケージ ‘AnnotationDbi’, ‘Biobase’, ‘BiocGenerics’, ‘BiocInstaller’, ‘BiocParallel’, ‘DelayedArray’, ‘GO.db’, ‘GenomeInfoDb’, ‘GenomeInfoDbData’, ‘GenomicRanges’, ‘HDF5Array’, ‘IRanges’, ‘Rhdf5lib’, ‘Rhtslib’, ‘S4Vectors’, ‘SingleCellExperiment’, ‘SummarizedExperiment’, ‘XVector’, ‘annotate’, ‘beachmat’, ‘biomaRt’, ‘bladderbatch’, ‘cellity’, ‘edgeR’, ‘genefilter’, ‘graph’, ‘impute’, ‘limma’, ‘org.Hs.eg.db’, ‘org.Mm.eg.db’, ‘preprocessCore’, ‘rhdf5’, ‘scPipe’, ‘scater’, ‘sva’, ‘topGO’, ‘tximport’, ‘updateR’, ‘zlibbioc’ が利用できません (for R version 3.5.0)
2:  install.packages(new, Ncpus = 6) で:
   一つもしくは複数のパッケージのインストールが失敗しました
       恐らく ‘WGCNA’

バージョンの確認

> version
               _
platform       x86_64-pc-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status
major          3
minor          5.0
year           2018
month          04
day            23
svn rev        74626
language       R
version.string R version 3.5.0 (2018-04-23)
nickname       Joy in Playing

インストールされているライブラリの確認

> nrow(installed.packages())
[1] 268

まぁいいでしょう。

参考:

https://makoto-shimizu.com/news/r-3-5-is-released/ https://askubuntu.com/questions/1031597/r-3-5-0-for-ubuntu https://launchpad.net/~marutter/+archive/ubuntu/rrutter3.5

googledrive内のファイルをRから操作 - googledrive

An Interface to Google Drive • googledriveを使ってみた

特徴

  • ほとんどの関数はdrive_で始まっているおかげでgoogledriveパッケージの関数を自動補完で呼び出しやすい。
  • find, ls, mv, cp, mkdir, rmといったUINIXコマンドを打つようにgoogleドライブを操作することを目標としている。
  • googledriveでは、"Drive tibble" 略してdribbleというオブジェクトにメタデータを保存する。dribbleはファイル名、ファイルIDなどを含む。
  • %>%でパイプが利用可能。

インストール

CRANからインストール可能

install.packages("googledrive")

使用法

まずは例の如くライブラリを読み込む。

library("googledrive")

ファイルの検索

drive_find(n_max = 50)
#> Auto-refreshing stale OAuth token.
#> # A tibble: 50 x 3
#>    name                    id                               drive_resource
#>  * <chr>                   <chr>                            <list>        
#>  1 chicken-xyz.csv         0B0Gh-SuuA2nTVUZGclZiSzZ0bkE     <list [37]>   
#>  2 chicken-rm.txt          0B0Gh-SuuA2nTT3dBbXd1ZWtvSkE     <list [38]>   
#>  3 chicken.jpg             0B0Gh-SuuA2nTbEhtYnIzcFNfX3M     <list [40]>   
#>  4 README-mirrors.csv      1LJlt-1emr662GV8WdEzddzsfqrt-Vg… <list [33]>   
#>  5 README-mirrors.csv      1PLXfempSnjpXbKVEXwMG5vBEnd-Fwm… <list [33]>   
#>  6 def                     0B0Gh-SuuA2nTRG5YWFVGaV8zbU0     <list [31]>   
#>  7 abc                     0B0Gh-SuuA2nTT2NqTGdLVWFkcjA     <list [31]>   
#>  8 folder1-level4          0B0Gh-SuuA2nTaTR6elE0TjZUUHM     <list [32]>   
#>  9 folder1-level3          0B0Gh-SuuA2nTWktWeTB0ajVoQjQ     <list [32]>   
#> 10 cranberry-TEST-drive-ls 1PM--xCb5axy5Uu9f6fDNjPAN2psRbQ… <list [32]>   
#> # ... with 40 more rows

パターンマッチ、ファイルタイプ(MIME)で検索をかけることもできる。

drive_find(pattern = "chicken")
drive_find(type = "spreadsheet")     ## Google Sheets!
drive_find(type = "csv")             ## MIME type = "text/csv"
drive_find(type = "application/pdf") ## MIME type = "application/pdf"

Google Drive API で利用可能なクエリパラメータを用いて検索することも可能。

'horsebean'という文字列をファイル内に含むものを検索したい場合は以下のようにする。

(files <- drive_find(q = "fullText contains 'horsebean'"))
#> # A tibble: 8 x 3
#>   name                             id                       drive_resource
#> * <chr>                            <chr>                    <list>        
#> 1 chickwts                         0B0Gh-SuuA2nTN05CNjk3bG… <list [38]>   
#> 2 chickwts_gdoc-TEST-drive-publish 1BHmmAyclG4RQS7hOJpeKlI… <list [32]>   
#> 3 foobar                           1qoA3kr9DmSTtsG9hoicP7y… <list [32]>   
#> 4 foobar                           0B0Gh-SuuA2nTa01CaXZZOW… <list [39]>   
#> 5 chickwts-TEST-drive-publish      0B0Gh-SuuA2nTSjh6SElzSV… <list [39]>   
#> 6 chickwts_gdoc-TEST-drive-list    1lA8iYCyFyFi7T_qe7UrYyB… <list [32]>   
#> 7 chickwts_txt-TEST-drive-publish  0B0Gh-SuuA2nTaU1CRmR2ZG… <list [38]>   
#> 8 hadley-googledrive-tour          1CndPuuAlTGNJkyqqE-CNuw… <list [32]>

ファイルのダウンロード

google drive中のファイルはGoogle Documents, Google Sheets, Google Slidesなどのnativeなファイルタイプになっているので、一般的なフォーマットに変換する必要がある。
デフォルトでファイルのフォーマットは自動的に選択されるが、明示的に指定する場合はtype=で指定する。

drive_download("538-star-wars-survey", 
                           type = "csv",
                           overwrite = TRUE)
#> File downloaded:
#>   * 538-star-wars-survey
#> Saved locally as:
#>   * 538-star-wars-survey.csv

path=を指定してファイルフォーマットを指定することも可能。

drive_download(
  "538-star-wars-survey",
  path = "538-star-wars-survey.csv",
  overwrite = TRUE
)
#> File downloaded:
#>   * 538-star-wars-survey
#> Saved locally as:
#>   * 538-star-wars-survey.csv

ファイルのアップロード

ファイルのアップロードにはdrive_upload関数を用いる。

> (chicken <- drive_upload(
  drive_example("chicken.csv"),
  "README-chicken.csv"
))
#> Local file:
#>   * /Users/jenny/resources/R/library/googledrive/extdata/chicken.csv
#> uploaded into Drive file:
#>   * README-chicken.csv: 1w9vv35y7pNE6q_wl3EYoh5ty_GO44gj0
#> with MIME type:
#>   * text/csv
#> # A tibble: 1 x 3
#>   name               id                                drive_resource
#> * <chr>              <chr>                             <list>        
#> 1 README-chicken.csv 1w9vv35y7pNE6q_wl3EYoh5ty_GO44gj0 <list [37]>

typeはデフォルトだと自動的に決められるが、明示的に指定することも可能。

> chicken_sheet <- drive_upload(
  drive_example("chicken.csv"),
  "README-chicken.csv",
  type = "spreadsheet"
)
Local file:
  * /home/kimoton/R/x86_64-pc-linux-gnu-library/3.4/googledrive/extdata/chicken.csv
uploaded into Drive file:
  * README-chicken.csv: 10IL80OCfOUrVNmQiyeXyz92ss_uko2zfDX5MxpAX-NA
with MIME type:
  * application/vnd.google-apps.spreadsheet

ファイルの公開

Google Documents, Google Sheets, Google Presentationのファイルは、drive_publish()関数で公開できる。

> (chicken_sheet <- drive_publish(chicken_sheet))
#> Files now published:
#>   * README-chicken.csv: 12ZzU5hS7GFQdpJBoVzBz3p8Xie2uNmCuNJNUVemGF9k
#> # A tibble: 1 x 7
#>   name     published shared id             drive_resource permissions_res…
#> * <chr>    <lgl>     <lgl>  <chr>          <list>         <list>          
#> 1 README-… TRUE      TRUE   12ZzU5hS7GFQd… <list [33]>    <list [2]>      
#> # ... with 1 more variable: revision_resource <list>

ファイルの公開状況を確認する際はdrive_reveal()を使う。

chicken_sheet %>% 
  drive_reveal("published")
#> # A tibble: 1 x 7
#>   name     published shared id             drive_resource permissions_res…
#> * <chr>    <lgl>     <lgl>  <chr>          <list>         <list>          
#> 1 README-… FALSE     TRUE   12ZzU5hS7GFQd… <list [33]>    <list [2]>      
#> # ... with 1 more variable: revision_resource <list>

ファイルの削除

drive_rmコマンドで行う。単純に引数に削除したいファイルを渡すだけ。

drive_rm(chicken_sheet, text_file)

参考

An Interface to Google Drive • googledrive

WSL (Windows Subsystem for Linux) で文字化け

An Interface to Google Drive • googledriveこれ使ってgoogledriveをコマンドからいじれるようにしたかったのにそもそもWSL内にgoogle-chromeが入ってなかった。

google-chromeをインストール

google-chromeをインストールします。

署名鍵のダウンロード、登録

sudo sh -c 'echo "deb http://dl.google.com/linux/chrome/deb/ stable main" >> /etc/apt/sources.list.d/google.list'
sudo wget -q -O - https://dl-ssl.google.com/linux/linux_signing_key.pub | sudo apt-key add -

これによりパッケージマネージャがGoogle Chromeの完全性を検証できるようになる。

install

今回はstable版をインストール

sudo apt-get install google-chrome-stable

起動

google-chrome

f:id:kimoppy126:20180608081827p:plain

もじばけ!

$ sudo apt-get install unifont

日本語フォントがインストールされてなかっただけでした。

f:id:kimoppy126:20180608082002p:plain

クラスタの特徴を知る - radarchart

クラスタリングを行ったあと、各クラスタがどんな特徴を持っているのか知りたいときはレーダーチャートを書いたりする。 radarchart関数を使う

radarchart(df)でとりあえずのレーダーチャートはかける。ここで与えるデータフレームは

df
The data frame to be used to draw radarchart. If maxmin is TRUE, this must include maximum values as row 1 and minimum values as row 2 for each variables, and actual data should be given as row 3 and lower rows. The number of columns (variables) must be more than 2.
radarchart function | R Documentation

1行目がmax、2行目がmin、3行目からが実際に可視化するデータとなっていなければならない。

kmeansを実行したデータフレームに適用する例

state.km <- kmeans(scale(state.x77[,1:6]),3)
df <- as.data.frame(scale(state.km$centers))
dfmax <- apply(df, 2, max)
dfmin <- apply(df, 2, min)
df <- rbind(dfmax, dfmin, df)

f:id:kimoppy126:20180527131910p:plain

ただこれだと、データのmax、minがレーダーチャートのmax、minとなってしまうため少し見づらい。

df.maxに +1、df.min に-1してやると、見やすくなる。

dfmax <- apply(df, 2, max) + 1
dfmin <- apply(df, 2, min) - 1
df <- rbind(dfmax, dfmin, df)

f:id:kimoppy126:20180527131917p:plain

まめち。

改訂2版 データサイエンティスト養成読本 [プロになるためのデータ分析力が身につく! ] (Software Design plus)

改訂2版 データサイエンティスト養成読本 [プロになるためのデータ分析力が身につく! ] (Software Design plus)

非線形クラスタリング k-means

k-meansの必要なところだけ。

k-meansとは

非階層的クラスタリング手法の1つ。

要するに何をしているのか

  1. k個のクラスターの初期位置を決める。
  2. 各データをk個のクラスターとの距離を求め、最も近い位置のクラスターに分類。
  3. 形成されたクラスターの中心を算出
  4. クラスタの中心が変化しない時点まで2, 3を繰り返す

参考:以下のデモツールをカチカチ https://www.albert2005.co.jp/knowledge/data_mining/cluster/non-hierarchical_clustering

短所

k-means法の一つの短所として、初期値(初期に選択される「核」となるk個のサンプル)依存性があります。図29の3つのクラスターは、初期値を変えて、重心が変化しなくなるまで、繰り返し計算した時の結果です。同じデータを距離などを同じ条件にして計算しても、初期値が異なるだけで、結果が大きく違うことが分かります。従って、よいクラスターを得るためには、初期値を変えて何回か分析を実施し、平均クラスター内距離が最小になる初期値を選択するなど、最適初期値での結果を採用することが望ましいといえます。

Rで実践

## PCA
state.pca <- prcomp(state.x77[,1:6], scale = T)

## kmeansでクラスタリング
state.km <- kmeans(scale(state.x77[,1:6]), 3)

## データフレームに変換、クラスターident列を追加。
state.pca.df <- data.frame(state.pca$x)
state.pca.df$cluster <- as.factor(state.km$cluster)

## Plot
ggplot(state.pca.df, aes(x = PC1, y = PC2, label = name, col = cluster)) +
    geom_text() +
    theme_bw()

スクリーンショット 2018-05-27 10.38.38.png (70.7 kB)

改訂2版 データサイエンティスト養成読本 [プロになるためのデータ分析力が身につく! ] (Software Design plus)

改訂2版 データサイエンティスト養成読本 [プロになるためのデータ分析力が身につく! ] (Software Design plus)

SAM format

SAM format

f:id:kimoppy126:20180413065958p:plain

リードをマッピングした結果の情報を示すのにSequence Alignment/Map (SAM) formatという形式がよく使われる。

SAM format は@から始まるheader行と、それに続く以下の11列から構成されるアライメントセクションによって構成される。 アライメントセクション中の各行はリード1本1本のマッピング結果を示している。

QNAME:
リードの名前

FLAG:
mapped/unmapped read, pairing,などを示すビットごとの数値。

RNAME:
リファレンス配列の名前。マッピングされなかったリードは * となる。

POS:
左から数えて最初にマッチした塩基のポジション。マッピングされなかったリードは0となる。

MAPQ:
マッピングクオリティー

CIGAR:
リファレンスに対し、どのようにアライメントされたかを示す。これは複数の要素を含む。各要素はoperator とnumberからなる。

MNAME:
paired-endの場合、もう一方のリードの名前を示す。"="の場合、もう一方のリードも同じ名前であることを示す。

MPOS:
最もよくマッピングされたポジションを示す。a number indicating the left most mapping position of the mate.

TLEN:
一つの要素からなる場合、他にも同様にマッピングされたリードが存在する場合、左側のリードの最もマップされた塩基から右側の最もマップされた塩基までの長さが計算される。

SEQ:
リファレンスとして使用されるリード中の配列。"*" でない場合、CCIGAR列のM/I/S/=/Xを足したものになる。

QUAL:
リードのクオリティー。fasta/q format中のものと同じ。

soft-clipping, hard-clipping

clipped の時点でアライメントはされたがくっつかなかったことを表す。その中で、くっつかなかった配列が、
→soft-clipping : リファレンス配列中に存在しない
→hard-clipping : リファレンス配列のどこかにその配列は存在する(キメラリードなど)。
いずれの場合もcoverageの計算には用いられるべきではない。

正しいかどうかは不明だけどこんな理解をしている。