なぜRをあえて「講義」で学ぶのか？

Rだけでバイオインフォマティクス解析を完結することは少ない。
シェルコマンドやPythonとの連携が一般的であり、全体の解析フローを構築する際にそれぞれのツールの役割を理解することが重要。

R vs Python: 二元論に意味はあるのか？

深層学習ではPythonが優位との意見も多いが、RとPythonは共に強力なコミュニティを持ち、それぞれの特性に強みと弱みがある。
二元論的な見方は、両者の本質を理解した場合には無意味。

AI時代のR学習: その意義とは？

AIツール（ChatGPTなど）を使えば、コード生成や解決策の提示が簡単に行える時代において、あえて講義形式で学ぶ意義は何か。
ネットで調べるか、AIツールに聞けば解決するケースが多い。

Rを学ぶ上で「講義形式」が重要な理由

初心者にとっては、ゼロからイチへの壁が大きい。
個別のコード学習では全体像が見えにくく、モチベーションを保てないケースも多い。
講義形式では、全体像の理解や具体的なイメージを持つことができる。

講義形式で学ぶことの意義: 体系的な知識と実践的スキル

講義形式での学びは、「体系的な知識の習得」と「実践的なスキルの習得」に焦点を当てる。
インターネットには断片的な情報が多いが、体系的に解析の流れを学び、データの取り扱いから解析まで一貫したスキルを身につけることが重要。

実際の解析シナリオでのRの活用法を学ぶ

Rのパッケージを活用し、データ操作や可視化を効率化。
tidyverse, ggplot2, dplyr, data.tableなどを組み合わせ、強力なデータ解析パイプラインを構築する。
パッケージの使い方を単独で学ぶだけでは、全体の流れを理解できない。

講義の目的

Rを通じてデータ解析の全体像を把握する
解析の全体像を理解し、Rを使った各ステップがどのように役立つかを学ぶ。
他のツール（シェルコマンドやPython）との連携を意識しながら、実践的な解析フローを理解する。

Rがバイオインフォマティクスで活用される場面

Rの役割と他ツールとの連携の重要性

バイオインフォマティクス解析において、Rはその強力な統計解析機能や可視化能力によって広範囲に活用される。
しかし、Rだけで全体の解析が完結することは稀であり、シェルコマンドやPythonなどのツールと連携しながら、解析フローを構築するのが一般的。
RNA-seq解析を例に、Rがどのように役立ち、他のツールとどのように連携するかを説明。

例：RNA-seq解析におけるRの具体的な役割

RNA-seq解析は、データの前処理から発現変動解析、結果の可視化まで、多段階にわたるフローが必要。
各ステップでRの果たす役割と、それを補完する他のツールを説明。

RNA-seqパイプラインにおける各処理ステップと対応ツール

処理ステップ	使用ツール	目的	スクリプト/プログラム言語
FastQファイルの結合	`cat`	リシーケンスされたFastQファイルをサンプルごとに結合	シェルスクリプト
FastQファイルのサブサンプリングとストランド推定	fq, Salmon	サブサンプリングしてライブラリのストランド情報を推定	Python, シェル
リードの品質管理 (QC)	FastQC	シーケンスリードの品質を評価	シェルスクリプト
UMI抽出	UMI-tools	UMIを抽出してPCR重複を補正	シェルスクリプト
アダプター除去と品質トリミング	Trim Galore!	アダプター除去と低品質リードのトリミング	シェルスクリプト
ゲノムコンタミナントの除去	BBSplit	ゲノムに由来するコンタミナントの除去	シェルスクリプト
リボソームRNAの除去	SortMeRNA	リボソームRNAの除去	シェルスクリプト
アライメントと定量化ルートの選択	STAR, Salmon, RSEM, HiSAT2	リードのアライメントと定量化	シェルスクリプト
アライメントのソートとインデックス作成	SAMtools	アライメントファイルのソートとインデックス作成	シェルスクリプト
UMIベースの重複除去	UMI-tools	UMIに基づいてPCRの重複リードを除去	シェルスクリプト
トランスクリプト組み立てと定量化	StringTie	トランスクリプトアセンブリと発現量の定量化	シェルスクリプト
広範な品質管理 (QC)	RSeQC, Qualimap, dupRadar, Preseq, DESeq2	発現データの品質評価と異常やバイアスの検出	R, シェルスクリプト
発現変動解析 \| DESeq2 \| 発現変動遺伝子の解析、統計モデルの構築 \| R \|
結果の可視化とQCのレポート	MultiQC, R	QC結果、発現量の可視化、サンプル間類似度の評価	R, シェルスクリプト

各ツールとスクリプトの役割

シェルスクリプト:
- 大量データの処理、ファイル操作、並列処理に使用され、ファイルの結合、トリミング、除去などの高速操作を担当。
Python:
- データの前処理、サブサンプリング、サンプルの準備段階で使用。柔軟なデータ操作と効率的なスクリプト実行が可能。
R:
- 最終的な解析と可視化において使用。特に発現変動解析（DESeq2）やQCの結果の可視化において重要な役割を果たす。

Rの具体的な役割

発現変動解析:
- ツール: DESeq2
- 発現変動解析の結果を導き出すために、Rで統計モデルを構築し、遺伝子発現の変動を評価。
結果の可視化:
- ツール: MultiQC, R
- QC結果、発現量の分布、サンプル間の類似度を可視化し、Rのグラフィックライブラリ（ggplot2など）を使用して直感的な図を作成。

Rがバイオインフォマティクスで活躍できる場面と限界のある場面

RNA-seq解析のフローはRが活躍する場面が限られているようにみえる
シェルやPythonで行う多くのタスクはRでも可能だが、以下のようにRのシステム的な制約や効率性を考慮すると、必ずしもRが最適な選択肢とは言えない
Rはシステム的な制約があり、必ずしも効率的でない場合もあるため、シェルやPythonと併用することが望ましい。

大規模データ処理におけるRの限界

Rは大規模データを扱う際、メモリ使用量や実行速度に制約がある。
数百ギガバイト単位のデータでは、シェルやPythonの方がパフォーマンスが優れている。
data.tableやdplyrなどを活用すれば、ある程度のパフォーマンス向上は可能だが、膨大なデータを効率的に操作するには、シェルやPythonの方が適している

Rの強みを最大限に活かす場面

Rは統計解析とデータ可視化に強みがあり、特にggplot2やshinyは他の言語にはない強み。
インタラクティブなデータ分析やカスタムメイドのデータ可視化が必要な場面ではRが非常に有効。

適材適所でツールを使い分ける

バイオインフォマティクス解析では、各ツールを特性に応じて使い分ける必要がある。
- シェル: ファイル操作や大規模データの処理、ソフトウェアに基づくパイプライン構築
- Python: 柔軟なデータ操作や、深層学習を中心とした機械学習
- R: 統計解析、機械学習や可視化

解析フローの設計

ツール間のデータ形式の整合性を保ち、統一した解析パイプラインを構築することが重要。
Rで前処理を行い、Pythonでモデル化するなど、ツールの強みを活かす。

R、シェル、Pythonを組み合わせる際のベストプラクティス

適材適所でツールを選択:
- データ前処理やファイル操作、パイプライン作成にはシェル、深層学習を中心とした機械学習、大規模データ処理にPython、統計解析や可視化にはRを使用。
再現性の確保:
- snakemakeやnextflowを使い、解析パイプラインの再現性を確保。
- rmarkdownやjupyter notebookで解析結果をドキュメント化。
環境管理:
- condaやDockerを使用して依存関係を管理し、解析環境の統一と再現性を確保。

Rの役割は下流解析において重要

生データを定量値に変換するプロセスでは、Rの役割は比較的小さい場合が多い。これは、多くの解析パイプラインで前処理や品質管理が自動化されているためである。
下流解析では、Rが重要な役割を果たす。複数のデータの統合、アノテーション、メタ情報の処理などを行い、定量されたデータを知識へと変換するプロセスで活躍。

Rを用いてできることの全貌

項目	説明	具体例
データの入出力と前処理	CSV、Excel、JSON、SQLなどのデータを効率的に読み込む。欠損値処理やスケーリングも容易。	`read.csv`, `fread`, `tidyverse`, `janitor`
データの操作と変換	データのフィルタリング、集計、結合やピボット操作が可能。文字列操作も効率的に行える。	`dplyr`, `data.table`, `tidyr`, `stringr`
データの可視化	基本的なプロットやカスタマイズ可能な高度な可視化が可能。インタラクティブなグラフ作成もサポート。	`plot`, `ggplot2`, `plotly`, `shiny`
統計解析とモデリング	基本的な統計解析や回帰モデル、ベイズ統計も扱える。混合効果モデルやサバイバル分析も可能。	`t.test`, `aov`, `glm`, `lme4`, `rstan`, `brms`
機械学習とAI	教師あり学習、教師なし学習、深層学習が可能。多くの機械学習アルゴリズムに対応。	`caret`, `randomForest`, `e1071`, `keras`, `tensorflow`
時系列解析	ARIMAモデルなどを使った時系列解析や異常検知が可能。	`forecast`, `ts`
空間データ解析	地理空間データの取り扱いや解析、地図の描画が可能。	`sf`, `sp`, `raster`
テキストマイニングとNLP	テキストデータの前処理やトピックモデリング、感情分析、文章分類などを行える。	`tm`, `text2vec`, `LDA`
リポート作成とプレゼン	`rmarkdown` を使った自動レポート生成や、`shiny` を用いたインタラクティブなダッシュボード作成が可能。	`rmarkdown`, `shiny`
プログラミングとアルゴリズム	基本的なプログラミング操作に加え、Rパッケージの作成や関数の共有が可能。	`devtools`, `usethis`

Rを用いた下流解析

下流解析には、標準的な解析を実行するためのRパッケージが多数存在。
Bioconductorは、R環境でのバイオインフォマティクス解析を支援するための大規模なソフトウェアプロジェクト
特にバイオインフォマティクス解析で利用され、次世代シーケンシング（NGS）やマイクロアレイデータの解析に特化。
各ステップで適切な関数を組み合わせて、全体の解析フローを構築する必要がある。

Bioconductorパッケージを使うとできること

分析タイプ	説明	主なパッケージ/ツール
遺伝子発現解析	マイクロアレイやRNA-seqデータの発現変動解析。パスウェイ解析や時系列解析も対応。	`affy`, `limma`, `DESeq2`, `edgeR`, `tximport`, `DEXSeq`, `clusterProfiler`, `ReactomePA`
エピゲノム解析	ChIP-seq、ATAC-seqデータのピークアノテーションや発現変動解析。 \| `ChIPseeker`, `DiffBind`, `csaw`, `chromVAR`, `ATACseqQC`
遺伝子構造と変異解析	ゲノムアノテーションとバリアント解析。構造変異や突然変異シグネチャ解析も対応。	`GenomicRanges`, `rtracklayer`, `biomaRt`, `VariantAnnotation`, `StructuralVariantAnnotation`, `MutationalPatterns`
メチル化解析	DNAメチル化データやバイサルファイトシーケンスデータの解析。	`minfi`, `RnBeads`, `BiSeq`, `methylKit`
マイクロバイオーム解析	メタゲノムや16S rRNAデータのインポート、正規化、可視化、解析。	`phyloseq`, `metagenomeSeq`, `dada2`, `microbiome`
シングルセル解析	シングルセルRNA-seqデータの前処理、クラスタリング、発現変動解析。	`scater`, `scran`, `SingleCellExperiment`, `Seurat`, `Scanpy`, `monocle`, `MAST`, `tradeSeq`
質量分析データ解析	メタボロームやプロテオームデータのピーク検出、統合解析、タンデム質量スペクトル解析。	`XCMS`, `MetaboAnalystR`, `CAMERA`, `MSnbase`, `ProteoMM`, `mzR`, `proFIA`, `Spectra`, `mzIdentML`
ネットワーク解析	遺伝子共発現ネットワークや相互作用ネットワークの解析と可視化。	`WGCNA`, `graph`, `STRINGdb`, `biogrid`, `pathview` \|
インタラクティブなデータ解析	Webアプリケーションやシングルセルデータのインタラクティブな可視化と操作。	`shiny`, `iSEE`
パイプラインとワークフロー	再現性の高い解析ワークフローの構築とデータ管理。	`BiocWorkflowTools`, `workflowr`, `Drake`, `ExperimentHub`, `AnnotationHub`, `MultiAssayExperiment`
知識ベースとアノテーション	ゲノムや遺伝子のアノテーション、アノテーションデータの取得と管理。	`biomaRt`, `GenomicFeatures`, `org.Hs.eg.db`, `AnnotationDbi`, `EnsDb`

Bioconductorの特筆すべき利点

データの一貫した取り扱い:
- SummarizedExperimentやSingleCellExperimentなどの統一されたデータ構造を使用することで、解析の各ステップ間でのデータ受け渡しがスムーズになる。
豊富なアノテーション情報の統合:
- biomaRtやAnnotationHubで簡単にアノテーション情報を取得し、解析結果の解釈を容易にする。
最新の解析手法への対応:
- RNA-seqやシングルセル解析に対応したツールを迅速に提供し、最先端の解析を行える。
コミュニティとサポート:
- 活発なコミュニティに支えられ、最新の技術やパッケージが常にアップデートされる。
再現性の高い解析:
- 並列処理やパッケージ管理機能により、解析結果の再現性を確保。

R + Bioconductorの強力な組み合わせ

RNA-seq解析やシングルセルRNA-seq解析では、データのインポート、解析、アノテーション付加、可視化を一貫して行える。
Rの可視化ツール（ggplot2やComplexHeatmap）を用いた高度な可視化が可能。
Bioconductorと連携することで、Rはバイオインフォマティクス解析の強力なツールとなる。

Pythonとの比較

Bioconductor相当のプラットフォームの不在
- Pythonには、RのBioconductorのような包括的なバイオインフォマティクス解析エコシステムがない。
統計解析とモデリングの強力なサポート
- Rは、高度な統計的手法（GLM、ベイズ統計など）の適用が容易で、RNA-seq解析では特に優れている。
バイオインフォマティクス特化のエコシステム
- Bioconductorは1000以上のパッケージで、RNA-seqやChIP-seqなど様々な分野に対応し、統合的な解析が可能。
データの標準化と再利用性の高さ
- SummarizedExperiment や SingleCellExperiment などのデータ構造が標準化され、パッケージ間でのデータのやり取りがスムーズ。
可視化の柔軟性と高品質なグラフ生成
- Rのggplot2は、ボルケーノプロットやヒートマップなどの複雑な可視化が容易で、高品質な図表が生成可能。
豊富なサポートとコミュニティ
- Bioconductorには活発なコミュニティがあり、ドキュメントやチュートリアルが充実している。

Rを用いた下流解析で必要な操作スキルとは？

RおよびBiocondutorパッケージを用いて下流解析フローを適切に設計して実装するためには、大まかに二つの種類の本質的理解が必要となる。

各パッケージ関数の背後にある理論的な理解
各パッケージ関数を活用する際に生じる前後のデータ整理と変換・可視化

本講義ではRの基礎を扱うため、1については扱わない。

簡単な補足：各関数がどのようなアルゴリズムや統計手法に基づいているかを理解することは、信頼性の高い解析を行う上で不可欠。
例えば：RNA-seqにおけるlimmaやDESeq2の発現変動解析では、正規化、分散のモデリング、仮説検定といった理論が関わる。
理論的背景を理解していないと、パラメータ設定や結果の解釈に誤りが生じる可能性がある。
体系的知識は、トピックごとに扱うことが具体的に適している。

主に2についてさらに考えていく。

Bioconductorのパッケージには、Vignetteという使用例がMarkdown形式で添付されていることが多い。
Vignetteに従うことで、解析フローを比較的簡単に実行することが可能。
初心者ユーザーでも、Vignetteに従うことで、関数に入力データを渡し、出力を次のステップに引き渡すといった基本的な操作が可能。
実際のデータはVignette通りの形式や条件で得られることは少ない。
例外的なケースに遭遇することが多く、その際、ユーザーはRの操作に習熟して対応する必要がある。

状況に対応するためのスキル

Rでパッケージ関数を組み合わせて解析フローを構築する際に必要なスキル：

関数を正しく使う
- 目的に応じた関数を選び、適切に使う技術が必要。
関数に対する入力データと引数を用意する
- 関数が要求する入力データやパラメータを整形・設定する必要がある。
関数からの出力データを整理・可視化する
- 出力を整形し、次のステップに入力したり、可視化・解釈できる形式に変換する。

Rでの解析における「親切な関数」と「不親切な関数」の存在

上記のスキルがなぜ必要か？

親切な関数
- 入力データがある程度整っていれば、自動的に多くの処理を行う。
不親切な関数
- データの整形や前処理をユーザーに要求する。
  - 開発者の立場からすると、親切な関数を作ることは非常に難しい。
  - あらゆる可能性に対応するのは膨大な作業量を必要とする。
  - ユーザーのスキルやニーズは多様で、どのように関数を使うかも予測しづらい。
  - ユーザーが自ら操作を行う必要が出てくるのは仕方がないことかもしれない。

関数の出力データの分析・整理の必要性

ユーザー自身がフィルター条件を設定し、出力データを分析する必要がある。
パッケージには出力データを可視化する関数が用意されているが、その多くは見栄えが悪いことがある。
カスタマイズの必要性から、別の可視化パッケージを使用する選択肢も出てくる。
この場合、ユーザーはその可視化のために出力データを整形・変換する必要がある。

必要なスキルを習得するためのステップ

以下に、RNA-seqの発現変動解析を例に取る。

各ステップでユーザーがどの操作を行うべきかを示す。
これにより、データ構造の理解やベクトル操作、データの整形・可視化の重要性が明確になる。

例1: RNA-seq発現変動解析（limma）

RNA-seqデータの発現変動解析のフローを示し、どの部分が自動的に処理され、どこでユーザーの介入が必要か解説。

# 1. データの読み込みとDGEListオブジェクトの作成
# ユーザーが行う部分
library(limma)
library(edgeR)

count_data <- read.table("counts_matrix.txt", header=TRUE, row.names=1)
meta_data <- read.csv("metadata.csv", row.names=1)

dge <- DGEList(counts=count_data, samples=meta_data)

# 2. アノテーションデータの統合
# ユーザーが行う部分
annotation_data <- read.table("gene_annotation.txt", header=TRUE, row.names=1)
dge$genes <- annotation_data[rownames(dge), ]

# 3. 正規化
# edgeRが自動的に行う部分
dge <- calcNormFactors(dge)

# 4. メタ情報を整理してデザインマトリクスを作成
# ユーザーが行う部分
group <- factor(meta_data$condition)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)

# 5. 遺伝子のフィルタリング
# ユーザーが行う部分
keep <- filterByExpr(dge, design)
dge <- dge[keep, , keep.lib.sizes=FALSE]

# 6. 分散縮小モデルの適用
# voomが行う部分
v <- voom(dge, design, plot=TRUE)

# 7. フィッティングと経験ベイズによる発現変動解析
# limmaが自動的に行う部分
fit <- lmFit(v, design)

# 8. コントラストの設定と解析実行
# ユーザーが行う部分
contrast <- makeContrasts(Treatment-Control, levels=design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

# 9. 結果の整形と出力
# ユーザーが行う部分
results <- topTable(fit2, adjust="fdr", number=Inf)
significant_genes <- results[results$adj.P.Val < 0.05 & abs(results$logFC) > 2, ]
write.csv(significant_genes, file="limma_differential_expression_results.csv")

各ステップの説明

DGEListオブジェクトの作成
- ユーザーがカウントデータとサンプル情報を読み込み、DGEListオブジェクトにまとめる。
アノテーションデータの統合
- ユーザーが外部のアノテーションデータを読み込み、DGEListオブジェクトに統合する。
正規化
- edgeRが自動的にサンプル間のライブラリサイズを補正する。
メタ情報の整理とデザインマトリクスの作成
- ユーザーがサンプル条件に基づき、デザインマトリクスを作成。
遺伝子のフィルタリング
- ユーザーが解析対象となる遺伝子をフィルタリングする。
分散縮小モデルの適用
- voom関数が自動的にモデルを適用。
フィッティングと経験ベイズによる発現変動解析
- limma関数が自動的に解析を行う。
コントラストの設定と解析実行
- ユーザーが比較条件を指定し、統計モデルを構築する。
結果の整形と出力
- ユーザーが結果を整形し、フィルタリング基準に基づいて重要な遺伝子を選別する。

定量アッセイに対する下流解析の一般的フローと必要なR操作

定量アッセイに対する下流解析の流れには、どのアッセイ系でも一定の「型」や「モデル」が存在する。
アッセイ系が多様でも、シークエンサーや質量分析（かつてはマイクロアレイ）の原理に依存しているため、共通する部分がある。
下流解析の本質を理解する上で、このメタ的な視点は重要。
測定限界やアーチファクト、ノイズ構造は多くの場合、異なるプロトコルであっても数理的には共通していることも多い。

1. データの統合とアノテーション

多くの定量アッセイでは、異なる実験やデータセットの統合が必要。
データセットの統合
- 異なるバッチやサンプルから取得したデータを統合。
- merge()やdplyrパッケージのleft_join()で共通のキー（遺伝子IDやタンパク質ID）に基づき統合。
アノテーションの追加
- 外部アノテーションデータ（HGNCやEnsemblデータベース）から遺伝子名や機能的アノテーションを取得。
- biomaRtやAnnotationDbiを使用。
メタ情報の処理
- サンプルや条件に関するメタ情報（性別、年齢、実験条件）を統合。
- SummarizedExperimentやSingleCellExperimentのcolData()関数を使って統合。

2. データの前処理

正規化やフィルタリング、バッチ効果の除去などは品質保持に重要。
データの正規化とスケーリング
- サンプル間のバッチ効果や測定バイアスを取り除く。
- DESeq2のrlog()やlimmaのvoom()でスケーリングと補正を行う。
フィルタリング
- サンプルの低品質データやノイズを除去。
- apply()関数で遺伝子ごとの分散を計算、フィルタリングを実施。
バッチ効果の補正
- 異なるバッチ間のバイアスを補正。
- ComBatやlimmaのremoveBatchEffect()を使用。

3. 次元削減とクラスタリング

次元削減
- 高次元データを2次元・3次元に圧縮し、パターンを視覚化。
- PCA、t-SNE、UMAPをSeuratやggplot2で可視化。
クラスタリング
- サンプルや条件をクラスタリングしてパターンを特定。
- K-meansや階層的クラスタリング、Seuratによるクラスタリング。

4. 機械学習と予測モデル

モデルの訓練とテスト
- 学習データで予測モデルを構築し、交差検証で性能を評価。
- caretやglmnetで複数のアルゴリズムを試し、最適なモデルを選定。
特徴選択
- 高次元データから重要な特徴を選択。
- 遺伝子発現データでは、発現変動遺伝子をモデルに使用。

5. 発現変動解析と機能的アノテーション

発現変動解析
- 条件間で遺伝子やタンパク質の発現量の違いを調べる。
- DESeq2やedgeRで統計的に有意な遺伝子を特定。
パスウェイ解析や機能的アノテーション
- 発現変動した遺伝子やタンパク質が関与するパスウェイや機能を明らかにする。
- clusterProfilerやReactomePAを使用。

6. ネットワーク解析

遺伝子やタンパク質の相互作用解析
- 相互作用ネットワークを構築し、重要な遺伝子やタンパク質を特定。
- WGCNAやigraphを使用してネットワーク解析を行う。

Rでデータを自在に操るための基礎

上記のデータ解析の「型」のうち、ユーザーが自ら行わなければならない本質的な操作に着目すると、データの前処理や整形、出力の可視化などの多くのステップが浮かび上がる。
以下は、そのような関数を利用する際に、ユーザーが習熟すべき主要なスキルをまとめたものである。

1. データの整形と前処理

データのインポートとフォーマットの調整
- データを適切な形式でインポートし、欠損値の処理や列のリネームなどを行うスキル。
- read.csv(), fread(), data.frame(), as.factor() などの基本操作。
データのクレンジング
- 欠損値や異常値を検出・修正し、データの質を確保するスキル。
- is.na(), na.omit(), dplyr::mutate() などを使用。

2. データのフィルタリングとサブセット化

データのサブセット化
- 必要なサンプルや変数だけを抽出し、解析用にデータをサブセット化するスキル。
- subset(), filter(), select() を使用。
フィルタリング
- 基準に基づいてデータをフィルタリングし、解析に適さない部分を除外するスキル。
- apply(), rowMeans(), filterByExpr() などを活用。

3. データのマージとアノテーション

データセットのマージ
- 異なるソースのデータセットを共通キー（遺伝子IDやサンプルID）で結合するスキル。
- merge(), dplyr::left_join() を使用。
アノテーションの追加
- 遺伝子やタンパク質のIDに関して、外部データベース（biomaRtやAnnotationDbi）からアノテーションを取得し統合する操作。

4. 関数の出力結果の整形と解析

出力データの整形
- 関数の出力データがリストや複雑な構造である場合に、必要な部分を抽出し、解析・可視化に使える形式に整形するスキル。
- unlist(), lapply(), sapply() などを使用。
結果のフィルタリング
- 統計的に有意な結果を抽出し、重要な結果を報告するスキル。
- topTable(), subset() でFDRやp値に基づくフィルタリングを実施。

5. カスタムの可視化

カスタムプロットの作成
- ggplot2やplotlyを使ってカスタマイズしたプロットを作成し、グラフの軸、凡例、タイトルなどを自由にカスタマイズするスキル。
結果の可視化
- クラスタリングや次元削減結果、発現変動解析の結果をプロットしてデータのパターンを視覚的に表現。
- PCA, t-SNE, UMAPを可視化するための技術を学ぶ (ggplot2, Seuratなどを使用)。

6. パイプラインの構築と自動化

データ処理のパイプライン化
- %>%（パイプオペレーター）を使った一連の操作を効率化。
- dplyrやtidyrを用いたパイプライン処理のスキルを習得。
スクリプトの自動化
- 同じ解析を繰り返す場合にスクリプトを作成し処理を自動化。
- forループやapply()を使って処理を繰り返すスキルを身につける。

パッケージ関数の活用を前提としたR操作の本質的な基礎とは？

上記の操作スキルには、R操作として重複している部分がある。それらをさらにまとめると、下流解析においてパッケージ関数の利用を前提としたRの基礎は、以下の点に集約される。

1. 入力・出力データのファイル形式およびデータ形式の多様性

生物学的データ解析では、さまざまなデータ形式やファイル形式（CSV、BAM、FASTQ、VCFなど）を取り扱う。
ユーザーは、Rのread.csv()やread.table()、RsamtoolsなどのBioconductorパッケージを使用してこれらのファイルを適切に読み込む必要がある。
データセットごとに異なるフォーマットを統一するため、データの変換や統合を行う操作スキルが必要。
ファイル形式の違いによって、読み込み時に適切な設定や前処理が必要であり、これらの操作はR操作の基本スキルに含まれる。

2. Rにおけるデータ構造の理解と操作スキル

Bioconductorパッケージでは、SummarizedExperiment、SingleCellExperiment、ExpressionSetといった特殊なS4クラスのデータ構造を扱う。
これらの構造に慣れるため、data.frameやmatrixに対する基本操作の理解が不可欠。
ユーザーは、基本的なデータフレーム操作に加え、S4オブジェクトのメタ情報（例：colData()やrowData()）を抽出する操作や、アッセイデータの取り扱いに習熟する必要がある。
これらのデータ構造に対する基本的な理解が、データ解析の前提となる。

3. データの整形とサマリー作成

Bioconductorパッケージの関数を使う前に、データを適切な形式に整形する必要がある。
遺伝子ごとの発現データやメタデータをフィルタリング・サブセット化し、解析に適した形に変換することが求められる。
発現変動解析を行う際には、入力データを正規化し、遺伝子やサンプルごとのフィルタリングが必要。
Rの基本的な関数（apply()やfilter()など）を使ってこれらの整形や変換操作を行う。
データの変換やサブセット化は、解析結果の品質に大きく影響するため、重要なスキルとなる。

4. 可視化のためのデータ整形と統計的要約

多くの解析結果は、最終的に可視化して解釈される。
主成分分析（PCA）やクラスタリング解析の結果を適切に可視化するために、データを整形・要約する必要がある。
ggplot2を使った可視化では、データをロング形式（long format）に変換し、必要な統計的要約を行うことが求められる。
適切な可視化を行うには、データ整形のスキルが欠かせない。また、プロットのカスタマイズやデータのサブセット化も重要なR操作の一部。

これらの4点は、Bioconductorパッケージ関数を活用する際の前提となる基礎的なR操作に該当。
パッケージ関数が提供する自動処理部分と、ユーザーが手動で行うデータ準備や変換・整形部分の間にある「操作の溝」を埋めるために、これらのRスキルが不可欠。
これらのスキルは解析の柔軟性を高め、問題が発生した際のトラブルシューティングにも役立つ。

データのライフサイクルを理解する

Rにおけるデータ操作を自在に行うためには、データのライフサイクルを把握することが基本となる。
バイオインフォマティクス解析において、データの流れは以下のサイクルに基づいて進行する。

データの読み込み
- CSV、Excel、FASTQ、BAMなど、多様な形式のデータをRに取り込む。
前処理
- データのクリーニング、欠損値の補完、データの正規化など、解析の前提条件を整える。
モデリングと分析
- 統計モデルや機械学習を用いて、データの構造を探り、予測や分類を行う。
要約と可視化
- 分析結果を要約し、視覚的にデータを理解するためのグラフやプロットを作成する。
レポート作成と結果の保存
- 解析の結果をまとめ、再現性のある形で記録・保存する。

このライフサイクルは一方向のプロセスではなく、前処理と分析、可視化を何度も繰り返しながら最適化を進めることが多い。

バイオインフォマティクス解析におけるデータの多様性

バイオインフォマティクス解析では、扱うデータの形式が非常に多様である。
一般的な表形式データだけでなく、配列データやゲノムアノテーションデータなど、独自の形式で保存されたデータを適切に読み込むための知識が不可欠。
次のようなデータ形式が存在する。

データ形式	説明	拡張子	Rパッケージ	読み込み関数	データ構造	書き込み関数	使用頻度
CSV/TSV	カンマまたはタブ区切りのテーブル形式データ	`sample.csv`, `sample.tsv`	`base R`, `data.table`	`read.csv()`, `fread()`	data.frame	`write.csv()`, `fwrite()`	高
Excel (XLS, XLSX)	Excelスプレッドシート	`sample.xlsx`, `sample.xls`	`readxl`, `openxlsx`	`read_excel()`, `read.xlsx()`	data.frame	`write.xlsx()`	中
JSON	JavaScript Object Notation。構造化データ	`sample.json`	`jsonlite`, `rjson`	`fromJSON()`, `toJSON()`	list/data.frame	`write_json()`	中
YAML	構造化データ（YAML形式）	`config.yaml`	`yaml`	`yaml.load()`, `yaml.load_file()`	list/data.frame	`write_yaml()`	低
XML	階層型データ	`data.xml`	`XML`, `xml2`	`xmlParse()`, `read_xml()`	list/data.frame	`write_xml()`	低
HDF5	階層型データフォーマット	`sample.h5`, `sample.hdf5`	`rhdf5`	`h5read()`, `h5ls()`	list/array	`h5write()`	中
RDS	Rオブジェクト保存形式	`object.rds`	`base R`	`readRDS()`	R object	`saveRDS()`	高
RData	複数のRオブジェクト保存形式	`data.RData`	`base R`	`load()`	environment	`save()`	中
Matrix Market (MTX)	スパース行列形式	`matrix.mtx`, `matrix.mtx.gz`	`Matrix`, `Matrix.utils`	`readMM()`	Matrix/Matrix-like	`writeMM()`	中
FASTQ	DNA/RNA配列リードシーケンスデータ	`sample.fastq`, `sample.fastq.gz`	`ShortRead`, `Biostrings`	`readFastq()`, `readDNAStringSet()`	ShortReadQ	`writeFastq()`	高
FASTA	配列データフォーマット	`sample.fasta`, `sample.fa`	`Biostrings`	`readDNAStringSet()`, `readAAStringSet()`	DNAStringSet	`writeXStringSet()`	高
SAM/BAM	シーケンスアライメントデータ	`sample.sam`, `sample.bam`	`Rsamtools`	`scanBam()`, `readBamFile()`	BamFile/BamFileList	`writeBamFile()`	高
VCF	バリアントコールフォーマット	`variants.vcf`, `variants.vcf.gz`	`VariantAnnotation`	`readVcf()`	VCF	`writeVcf()`	中
GFF/GTF	ゲノムアノテーションフォーマット	`annotation.gff`, `annotation.gtf`	`rtracklayer`	`import()`	GRanges	`export()`	中
BED	ゲノム領域データ	`regions.bed`	`rtracklayer`, `GenomicRanges`	`import()`	GRanges	`export()`	中
Count Matrix	遺伝子発現量のカウントデータ	`counts.txt`, `counts_matrix.tsv`	`edgeR`, `DESeq2`	`read.table()`, `read.delim()`	matrix/data.frame	`write.table()`, `write.delim()`	高
mzML/mzXML	質量分析のデータ	`sample.mzML`, `sample.mzXML`	`mzR`, `MSnbase`	`openMSfile()`, `readMSData()`	MSnSet	`writeMSData()`	中
BigWig	ゲノムカバレッジデータ	`coverage.bw`	`rtracklayer`	`import.bw()`	GRanges	`export.bw()`	中
AnnData (H5AD)	シングルセルデータ用HDF5形式	`sample.h5ad`	`zellkonverter`	`readH5AD()`	SingleCellExperiment	`writeH5AD()`	中
SQLite	SQLデータベース形式	`database.sqlite`	`DBI`, `RSQLite`	`dbConnect()`, `dbReadTable()`	SQLiteConnection	`dbWriteTable()`	低
Feather	高速読み込み用バイナリ形式	`data.feather`	`feather`, `arrow`	`read_feather()`, `write_feather()`	data.frame	`write_feather()`	低
Parquet	高速読み込み用列指向フォーマット	`data.parquet`	`arrow`	`read_parquet()`, `write_parquet()`	data.frame	`write_parquet()`	中
GZ/BZ2	圧縮形式（gzip/bzip2）で保存されたファイル	`data.gz`, `data.bz2`	`R.utils`, `base R`	`gunzip()`, `gzfile()`, `bunzip2()`	CompressedFile/Raw	`gzip()`, `bzfile()`	高
RLE	ランレングスエンコーディング形式	`data.rle`	`Rle`	`runLength()`	Rle	N/A	低
GMT	遺伝子セット定義フォーマット	`geneset.gmt`	`GSEABase`, `fgsea`	`getGmt()`, `gmtPathways()`	GeneSetCollection	N/A	中
SIF	ネットワーク相互作用データ形式	`network.sif`	`igraph`	`read.graph()`	igraph/network	`write.graph()`	低
TOML	構造化設定ファイルフォーマット	`config.toml`	`toml`	`toml::from_toml()`	list/data.frame	`toml::to_toml()`	低

データ読み込みの基盤スキル

データ読み込みのスキルの本質は、想定するデータ形式に沿わない場合に、適切な引数や関数のオプションを活用し、データを正確に読み込むための試行錯誤ができるスキルにある。
バイオインフォマティクスにおけるデータは非常に多様であり、必ずしも整った形式で提供されるとは限らない。そのため、以下のような能力が特に重要である。

1. データ形式の理解

さまざまなファイル形式（CSV, TSV, FASTQ, VCF, GTF, GMT, SIFなど）に対応し、それぞれの形式の特徴を理解することが第一歩である。
列の区切り文字、ヘッダーの有無、改行コードの違いなど、データ形式の違いが読み込みに与える影響を理解する必要がある。

2. 引数やオプションの柔軟な使用

各読み込み関数（fread, read.table, read.csv, read.delim など）には、データのフォーマットに合わせて引数を調整する機能が豊富に用意されている。
これらを活用して、データ形式に応じた適切な読み込みを行うことがスキルの中心である。
区切り文字の指定 (sep, delimiter): 区切り文字が予期せず異なる場合（例: カンマ区切り vs. タブ区切り）、引数を使って明示的に指定する。
```
fread("data.txt", sep = "\t")  # タブ区切りの場合
```
ヘッダーの有無 (header): ファイルの先頭にヘッダーがない場合や、複数行のヘッダーがある場合、header = FALSE や skip 引数を使う。
```
fread("data.txt", header = FALSE, skip = 2)  # ヘッダーが2行分ある場合
```
データ型の強制変換 (colClasses, as.is): 特定の列が予期せぬデータ型（例: 数値が文字列として扱われる）として読み込まれた場合、colClasses 引数を用いて強制的に型を指定する。
```
fread("data.txt", colClasses = c("GeneID" = "character", "Expression" = "numeric"))
```
欠損値の処理 (na.strings): 欠損値が異なる表現（例: “NA”, “-”, “null” など）で表記されている場合、それらをna.stringsで明示して読み込み時に処理する。
```
fread("data.txt", na.strings = c("NA", "-", "null"))
```

3. 試行錯誤とデータの確認

読み込み後のデータが期待通りであるかどうかを確認し、必要に応じて再調整するプロセスが重要である。
データの一部を確認するために、head() や str() などの関数を用い、誤って読み込まれた箇所があれば再度引数を調整して試行錯誤する。
```
# データの構造を確認
str(quant_data)
head(quant_data)
```

4. 想定外のデータフォーマットへの対応力

時には、データが標準的な形式ではなく、不規則な改行や複数の区切り文字が混在している場合もある。
このような場合、読み込み関数だけで対処できないことも多く、テキスト処理やファイルを前処理する技術が求められる。
正規表現を用いた前処理: テキストファイルが一貫性を欠いている場合、sed や awk などのコマンドラインツールや、Rの stringr パッケージで前処理を行う。
```
library(stringr)
data <- readLines("data.txt")
data_clean <- str_replace_all(data, "[^[:alnum:]\t]", "")  # 不要な記号を削除
```

5. データがどうしても読み込めない場合の対応

データが標準的な関数ではどうしても読み込めない場合、データを一旦 readLines や scan で読み込み、手作業でクリーニングすることがある。

このような手法を使うことで、複雑なデータ構造に対応できる。

# readLinesでファイルを行単位で読み込む
raw_data <- readLines("data.txt")

# 不要な行や空行を削除
raw_data_clean <- raw_data[raw_data != "" & !grepl("^#", raw_data)]  # 空行やコメント行を削除

# データを適切な形式に変換
data_matrix <- read.delim(textConnection(raw_data_clean), sep = "\t")

6. シェルスクリプトによる事前処理

Rだけでなく、シェルスクリプトを用いたデータ整形も有効である。
特に大規模なデータセットに対しては、シェルで事前に処理することで効率化できる。
シェルを使った基本的な整形操作を以下に示す。

1. `head` でデータの確認

先頭行からデータの形式を推測

head data.txt  # ファイルの先頭10行を確認

2. `sed` や `awk` を用いたデータ整形

区切り文字の変換（例: カンマをタブに置換）
```
sed 's/,/\t/g' data.txt > cleaned_data.txt
```

不要な行の削除（例: 空行や特定のパターンを含む行を削除）

sed '/^$/d' data.txt > cleaned_data.txt  # 空行を削除
awk '!/^#/' data.txt > cleaned_data.txt  # コメント行を削除

列の抜き出し（例: 2列目と4列目を抽出）
```
cut -f2,4 data.txt > selected_columns.txt
```

3. `grep` で特定行を抽出

特定の文字列（例えば、以下の”Pattern”という文字）を含む行を抽出
エラーメッセージなどから該当行を特定したり、除外（その行以外を抽出）したりする場合
特定の文字列（ヘッダーなど）を含む行のみ抽出
```
grep "Pattern" data.txt > filtered_data.txt
```

7. データ品質チェックとクリーニング

最後に、読み込んだデータが解析に適しているかどうかの品質チェックも不可欠である。
読み込み時のエラーや警告を無視せず、データを正しくクリーニングするプロセスを踏むことで、後続の解析をスムーズに進められる。
データ読み込みのスキルは、単に関数を使うだけでなく、データ形式の多様性に柔軟に対応し、問題に応じて引数を調整して試行錯誤する力にある。
特に、バイオインフォマティクスにおいては、データ形式の多様さと複雑さに応じて適切な処理を行い、シェルスクリプトなども活用して柔軟に対応することが求められる。
これらのスキルは、データ解析の最初のステップとして極めて重要であり、正確なデータの読み込みが解析結果の信頼性を高める基盤となる。

データ読み込みスキルの演習例題

以下に、データ読み込みのスキルを養うための例題をいくつか示す。これらの演習は、異なる形式のデータを読み込み、適切に整形するための技術を実践的に学ぶことが目的である。

例題1: CSVファイルの読み込みと整形

課題:

次のようなCSVファイル gene_expression.csv を読み込んでみよう。
このデータには、ヘッダー行が2行存在し、カンマ区切りで遺伝子発現量データが格納されている。
1行目はファイルの説明、2行目がヘッダーである。空白セルはNAとして扱う。

gene_expression.csv:

# Gene expression data from experiment XYZ
GeneID,SampleA,SampleB,SampleC
ENSG000001,10.5,9.8,,
ENSG000002,7.5,NA,12.3
ENSG000003,5.0,4.8,5.1

Rでのデータ生成と保存

# CSVファイルを生成し保存
gene_expression_data <- "GeneID,SampleA,SampleB,SampleC
ENSG000001,10.5,9.8,,
ENSG000002,7.5,NA,12.3
ENSG000003,5.0,4.8,5.1"
writeLines(c("# Gene expression data from experiment XYZ", gene_expression_data), "gene_expression.csv")

演習:

ファイルを読み込み、ヘッダー行をスキップして2行目からデータを扱うように設定する。
欠損値（空白セル）をNAとして処理する。
GeneID列を文字列型に指定し、数値列の型を確認する。
データの一部を表示し、正しく読み込めたかを確認する。

解答例:

# CSVファイルの読み込み
gene_data <- fread("gene_expression.csv", skip = 1, na.strings = c("", "NA"), fill = TRUE)
# GeneID列を文字列型にする
gene_data[, GeneID := as.character(GeneID)]

# データの確認
str(gene_data)
head(gene_data)

解説：

ステップ 1: CSVファイルの読み込み

gene_expression.csvファイルには2行のヘッダー行が含まれており、1行目はファイルの説明、2行目が実際のデータの列名を含んでいる。
したがって、fread()関数のskip = 1オプションを使用して1行目をスキップし、2行目からデータを読みこむ。
また、空白セル（欠損値）はNAとして扱うため、na.stringsオプションを設定

# CSVファイルの読み込み
gene_data <- fread("gene_expression.csv", skip = 1, na.strings = c("", "NA"), fill = TRUE)

skip = 1: 1行目をスキップし、2行目から読み込む設定。
na.strings = c("", "NA"): 空白セルやNAを欠損値として扱う設定。
fill = TRUE: 列の数が揃わない行がある場合に、自動的に空欄を埋める設定。

ステップ 2: `GeneID`列を文字列型にする

GeneID列は遺伝子IDを表すため、数値型ではなく文字列型として扱う必要がある
data.tableの列に対するデータ型の変更は、:=を使って行う

# GeneID列を文字列型にする
gene_data[, GeneID := as.character(GeneID)]

ステップ 3: データの確認

str()関数を使って、データフレームの構造を確認し、正しく読み込めたかをチェック
また、head()関数でデータの最初の数行を表示して、各列の内容が期待通りかを確認します。

# データの確認
str(gene_data)
head(gene_data)

例題2: タブ区切りファイルの読み込みと前処理

課題:

次のようなタブ区切りファイル sample_metadata.txt には、実験に関するメタデータが格納されている。
このデータにはいくつかの行がコメント（#で始まる行）として記載されており、空白行も含まれている。
このファイルを読み込み、不要なコメント行や空白行を削除してみよう。

sample_metadata.txt:

# Sample metadata from RNA-seq experiment
SampleID   Condition   Replicate
# Comment about the experiment
SampleA    Control     1
SampleB    Control     2

SampleC    Treatment   1
SampleD    Treatment   2

Rでのデータ生成と保存

# タブ区切りファイルを生成し保存
sample_metadata <- "# Sample metadata from RNA-seq experiment
SampleID\tCondition\tReplicate
# Comment about the experiment
SampleA\tControl\t1
SampleB\tControl\t2

SampleC\tTreatment\t1
SampleD\tTreatment\t2"
writeLines(sample_metadata, "sample_metadata.txt")

演習:

コメント行（#で始まる行）と空白行を削除し、データを読み込む。
各列のデータ型を確認し、正しく読み込めているか確認する。

解答例:

# ファイルを行単位で読み込む
raw_metadata <- readLines("sample_metadata.txt")

# 空白行とコメント行を削除
clean_metadata <- raw_metadata[raw_metadata != "" & !grepl("^#", raw_metadata)]

# データをタブ区切りで読み込む
metadata <- read.delim(textConnection(clean_metadata), sep = "\t")

# データの確認
str(metadata)

解説：

ステップ 1: ファイルを行単位で読み込む

readLines()関数を使ってファイルを行ごとに読み込みます。この関数はファイル内の各行をベクトルとして読み込むため、後で不要な行（コメント行や空白行）をフィルタリングしやすくなります。

# ファイルを行単位で読み込む
raw_metadata <- readLines("sample_metadata.txt")

ステップ 2: コメント行と空白行を削除する

コメント行は#で始まる行であり、grepl()関数を使ってそれらの行を検出する。
また、空白行も削除するために、行が空でないかを確認する。
このフィルタリングは、!= ""で空白行を除き、!grepl("^#", raw_metadata)で#で始まる行を除去することで行う

# 空白行とコメント行を削除
clean_metadata <- raw_metadata[raw_metadata != "" & !grepl("^#", raw_metadata)]

ステップ 3: データをタブ区切りで読み込む

clean_metadataには、不要な行が削除されたデータが含まれている
このクリーンなデータをタブ区切りで読み込むために、read.delim()関数を使う
textConnection()関数を使うことで、文字列データをファイルのように扱うことができる

# データをタブ区切りで読み込む
metadata <- read.delim(textConnection(clean_metadata), sep = "\t")

ステップ 4: データの確認

最後に、str()関数を使って、データフレームの構造を確認し、正しくデータが読み込まれているかをチェック。

# データの確認
str(metadata)

例題3: FASTQファイルの読み込みとクリーニング

課題:

次のようなFASTQ形式のファイル reads.fastq には、シークエンスリードデータが含まれている。
このファイルには改行が不規則な部分が含まれているため、まずはテキストとして読み込み、クリーニングする必要がある。
最終的に、各リードが4行ごとに正しく整列されているかを確認する。

reads.fastq:

@SEQ_ID1
GATTTGGGG
+
!!!!!
@SEQ_ID2
GATCAAGGA
+
!!!!!
@
SEQ_ID3
GATTAGGAA
+
!!!!!

リードは４行ごとに

1行目: リードのID（@SEQ_ID）
2行目: シークエンスデータ
3行目: “+” シンボル（品質スコアの開始）
4行目: 品質スコア（同じ長さのデータ）

となっている必要。しかし、ID3のみIDに破損があり、正しく読み込めない。

対応策：

空行や不適切な改行の削除
- @SEQ_ID3 の部分で見られる不規則な改行を修正する必要。
- 不規則な改行を削除し、各リードが4行セットに揃うようにクリーニングする必要がある。

Rでのデータ生成と保存

# FASTQ形式ファイルを生成し保存
fastq_data <- "@SEQ_ID1
GATTTGGGG
+
!!!!!
@SEQ_ID2
GATCAAGGA
+
!!!!!
@
SEQ_ID3
GATTAGGAA
+
!!!!!"
writeLines(fastq_data, "reads.fastq")

演習:

readLines()を使ってファイルを読み込み、リードの欠けている行を手作業で修正する。
各リードが4行ごとに正しい形式で並んでいるかを確認する。

解答例:

# FASTQファイルを読み込む
fastq_data <- readLines("reads.fastq")

# 空行を削除
fastq_data_clean <- fastq_data[fastq_data != ""]

# '@'で始まらない行を修正
for (i in seq_along(fastq_data_clean)) {
  if (grepl("^@", fastq_data_clean[i]) && i < length(fastq_data_clean) && !grepl("^@", fastq_data_clean[i + 1])) {
    fastq_data_clean[i] <- paste0(fastq_data_clean[i], fastq_data_clean[i + 1])
    fastq_data_clean[i + 1] <- ""  # 結合した次の行を空にする
  }
}
fastq_data_clean <- fastq_data_clean[fastq_data_clean != ""]  # 空行を再度削除

# リードが4行ごとに整列しているか確認
n_lines <- length(fastq_data_clean)
if (n_lines %% 4 == 0) {
  print("リードが正しく整列しています")
} else {
  print("リードの整列に問題があります")
}

# データを表示
print(fastq_data_clean)

解説：

ループの設定 (for (i in seq_along(fastq_data_clean))):
- seq_along(fastq_data_clean) は、fastq_data_clean ベクトルのインデックス番号を取得。つまり、fastq_data_clean に含まれる各行に対して順番に操作を行う。
条件判定 (if):
```
if (grepl("^@", fastq_data_clean[i]) && i < length(fastq_data_clean) && !grepl("^@", fastq_data_clean[i + 1]))
```
この条件は、以下の3つの条件を同時に確認している。
- grepl("^@", fastq_data_clean[i]):
  - i 番目の行が @ で始まるかどうかを確認。grepl("^@", ...) は正規表現を使って、文字列の先頭が @ で始まっているかをチェック。^ は「行頭」の意味を持つため、@ で始まる行を探している。
  - つまり、ここではシーケンスIDの行（@SEQ_ID）かどうかを確認している。
- i < length(fastq_data_clean):
  - i が fastq_data_clean の長さ（行数）を超えないことを確認している。これは、i+1 行目を参照する際に、範囲外アクセスを防ぐため。
- !grepl("^@", fastq_data_clean[i + 1]):
  - i+1 番目の行が @ で始まっていないことを確認。これは、シーケンスID行 (@SEQ_ID) が連続していないことを意味している。通常、シーケンスID行はリードの先頭で、次の行にはシーケンスが続くべきなので、@ が連続していないかを確認している。
行の結合 (paste0):
```
fastq_data_clean[i] <- paste0(fastq_data_clean[i], fastq_data_clean[i + 1])
```
- 上記の条件を満たす場合（つまり、@SEQ_ID の次に @ で始まらない行が続いている場合）、現在の行と次の行を結合する。paste0() は2つの文字列を結合する関数。
- 例えば、@SEQ_ID とその後のシーケンスが分かれていた場合に、それらを1行にまとめて正しいフォーマットに修正する。
次の行を空にする:
```
fastq_data_clean[i + 1] <- ""
```
- 結合した後、次の行 (i+1 行目) は不要になるので、空文字に置き換える。
空行の再削除:
```
fastq_data_clean <- fastq_data_clean[fastq_data_clean != ""]
```
- 空文字になった行をすべて削除して、リストをクリーンアップする。

例題4: CSVファイルの事前チェックとシェルスクリプトによる整形（参考）

課題:

次のような（仮にRで対処できない）大規模なCSVファイル large_data.csv を読み込む必要がある。
しかし、ファイルが非常に大きいため、Rに読み込まずに、まずシェルスクリプトを使ってファイルの先頭を確認し、その後、効率的なシェルコマンドで不要な列を事前に削除したいと考えている。
以下のステップに従って処理を進めてみよう。

large_data.csv:

SampleID,GeneA,GeneB,GeneC,GeneD,Metadata
1,5.4,3.2,7.1,2.4,SampleA_Info
2,6.1,2.8,6.9,3.3,SampleB_Info
3,5.9,3.5,7.4,2.9,SampleC_Info
...

データ生成と保存

# 大規模なCSVファイルを生成し保存
large_data <- "SampleID,GeneA,GeneB,GeneC,GeneD,Metadata
1,5.4,3.2,7.1,2.4,SampleA_Info
2,6.1,2.8,6.9,3.3,SampleB_Info
3,5.9,3.5,7.4,2.9,SampleC_Info"
writeLines(large_data, "large_data.csv")

演習:

シェルコマンドを使ってファイルの先頭を確認し、不要な列（Metadata列）を削除する。
整形されたファイルをRで読み込み、データの一部を確認する。

解答例:

シェルでの整形操作（Linux環境の場合）

# ファイルの先頭10行を確認
head large_data.csv

# Metadata列を削除して新しいファイルに保存
cut -d, -f1-5 large_data.csv > cleaned_large_data.csv

解説：

cut コマンドを使って、CSVファイルから不要な列を事前に削除することができる。
-d, は、CSVファイルのフィールド区切りがカンマであることを示す。
-f1-5 は、1列目から5列目まで（つまり、Metadata列を除いた部分）を選択して、新しいファイル cleaned_large_data.csv に出力する指示。

Rでの読み込み

# 整形されたCSVファイルを読み込む
large_data <- fread("cleaned_large_data.csv")

# データの確認
str(large_data)
head(large_data)

オブジェクトの構造を理解する重要性

Rにデータを読み込むと、そのデータはオブジェクトとして保存される（表中の「データ構造」を参照）。
オブジェクトの構造を正しく理解しないと、適切なデータ操作や解析が行えない。
- Rの関数は、特定のオブジェクト構造（ベクトル、行列、データフレーム、リストなど）を前提として動作する。
- 入力値や出力値の構造に誤りがあるとエラーが発生する。
特に、バイオインフォマティクス解析では、SummarizedExperiment、GRanges、SingleCellExperiment などの高度なデータ構造が頻繁に使われる。
- これらのオブジェクトを適切に操作するためには、構造を理解し、各ステップに応じた操作を選択する能力が必要。

オブジェクト間の変換と柔軟なデータ操作

解析を進める上で、データ構造を別の形式に変換する必要が生じることがある。
- 例: 行列形式で保持しつつ、解析結果をデータフレームに変換して可視化を行う。
Rには、as.data.frame() や as.matrix() などの関数があり、データの形式を適切に変換できる。
- これにより、異なるデータソースの統合や解析結果の変換が容易になる。
- 解析の効率が向上する。

継続的な試行錯誤とデータ解析の最適化

データのライフサイクルとオブジェクトの構造を理解した上で、Rによるデータ操作は一層自在になる。
探索的データ解析（EDA）では、前処理やフィルタリング、可視化を繰り返しながらデータの本質に迫る。
- 最適な解析手法を見つけ出す。
この過程では、Rの柔軟性を活かしつつ、試行錯誤を続けて解析の最適化を図ることが重要。
str() などの関数を使って、オブジェクトの詳細な構造を確認すると効果的。
- これにより、複雑なデータ構造を理解し、効率的にデータを操作できるようになる。

データ構造	種別	用途	データ構造を扱うパッケージ	サブセットの方法	データの基本処理のための主要な関数	適用可能な解析関数
ベクトル（vector）	CRAN	単一データ型の要素の集まり、基本的な統計計算やフィルタリング。	`base`	- `[` や `[[]]` でのインデックス指定 - 論理ベクトルでの条件指定	- 基本関数: `mean()`, `sum()`, `length()`, `sort()`, `unique()` - 演算関数: `+`, `-`, `*`, `/`, `^`	- 統計関数: `cor()`, `sd()`, `var()`, `t.test()`
リスト（list）	CRAN	異なるデータ型や構造をまとめ、複雑なデータの管理に使用。	`base`, `purrr`	- `[[ ]]` でのインデックス指定 - `$` での名前指定	- リスト操作: `lapply()`, `sapply()` (`base`) - リスト操作: `map()`, `map2()`, `pmap()` (`purrr`)	- リストの要素に関数を適用: `lapply()`, `map()`
行列（matrix）	CRAN	数値データの統計解析や行列演算。	`base`, `Matrix`	- `[ , ]` でのインデックス指定 - `row()` や `col()` で行・列の抽出	- 行列演算: `matrix()`, `%*%`, `t()`, `solve()` - 統計処理: `prcomp()` (`stats`)	- 主成分分析: `prcomp()` (`stats`) - 共分散行列: `cov()`
データフレーム（data.frame）	CRAN	異なるデータ型を持つ列を行列形式で管理し、データ解析の基本。	`base`, `dplyr`, `tidyr`	- `[ , ]` でのインデックス指定 - `subset()`, `filter()` (`dplyr`)	- データ操作: `mutate()`, `select()`, `filter()` (`dplyr`) - データ整形: `pivot_longer()`, `pivot_wider()` (`tidyr`)	- 統計解析: `lm()`, `glm()`, `aov()` (`stats`) - データ整形: `gather()`, `spread()` (`tidyr`)
ティブル（tibble）	CRAN	データフレームの改良版、`tidyverse` と連携しデータ解析に使用。	`tibble`, `dplyr`, `tidyverse`	- `[ , ]` でのインデックス指定 -`filter()`, `select()` (`dplyr`)	- データ操作: `mutate()`, `select()`, `filter()` (`dplyr`) - ティブル変換: `as_tibble()`, `glimpse()` (`tibble`)	- データ整形: `pivot_longer()`, `pivot_wider()` (`tidyr`) - データ解析: `summarize()`, `group_by()` (`dplyr`)
igraph	CRAN	グラフデータ構造（ネットワーク）の管理と解析。	`igraph`	- グラフのサブセット: `induced_subgraph()`, `delete_vertices()`	- ノード操作: `add_vertices()`, `delete_vertices()` - エッジ操作: `add_edges()`, `delete_edges()`	- ネットワーク解析: `betweenness()`, `closeness()`, `page_rank()` (`igraph`)
dgeMatrix（ディジーマトリックス）	CRAN	スパース行列の管理、行列演算や大規模データの解析に使用。	`Matrix`, `BiocParallel`	- `[ , ]` でのインデックス指定 - 行・列の名前指定: `rownames()`, `colnames()`	- 行列演算: `t()`, `%*%`, `solve()`, `crossprod()` - スパース行列の作成: `Matrix()`, `bandSparse()` (`Matrix`)	- 大規模行列の解析: `svd()`, `eigen()` - 主成分分析: `prcomp()`, `pca()`
SummarizedExperiment	Bioconductor	遺伝子発現データやアノテーション情報の統合管理。	`SummarizedExperiment`, `DESeq2`	- `assay()`, `rowData()`, `colData()` を用いた抽出	- データ操作: `assay()`, `rowRanges()`, `colData()` - フィルタリング: `rowRanges() <-`	- 発現変動解析: `DESeq()`, `results()` (`DESeq2`) \| - データ変換: `rlog()`, `vst()` (`DESeq2`)
ExpressionSet	Bioconductor	マイクロアレイやRNA-seqの発現データの統合管理。	`Biobase`, `limma`	- `exprs()`, `pData()`, `fData()` を用いた抽出	- データ操作: `exprs()`, `pData()`, `fData()` (`Biobase`) - データ整形: `addNA()`, `subset()`	- 発現変動解析: `lmFit()`, `eBayes()` (`limma`) \| - 発現変動解析: `voom()` (`limma`)
GRanges（GenomicRanges）	Bioconductor	染色体位置情報の管理、ゲノムアノテーションやピークデータの操作。	`GenomicRanges`, `rtracklayer`	- `GRanges()` オブジェクトのメタ情報を用いた抽出	- ゲノム領域の操作: `findOverlaps()`, `flank()`, `reduce()` - 領域間の操作: `subsetByOverlaps()`, `resize()`	- アノテーション解析: `annotatePeak()` (`ChIPseeker`) - 重複解析: `findOverlaps()` (`GenomicRanges`)
SingleCellExperiment	Bioconductor	シングルセルRNA-seqデータの統合管理。	`SingleCellExperiment`, `scran`	- `counts()`, `colData()`, `rowData()` を用いた抽出	- データ操作: `counts()`, `colData()`, `reducedDim()` - 正規化: `normalize()`, `logNormCounts()` (`scran`)	- 発現変動解析: `DEsingle()`, `MAST()` \| - クラスタリング解析: `quickCluster()`, `computeSumFactors()` (`scran`)
DGEList	Bioconductor	RNA-seqカウントデータの管理、発現変動解析の準備。 \| `edgeR` \| - `counts()`, `samples()`, `genes()` を用いた抽出 \| - データ操作: `calcNormFactors()`, `estimateCommonDisp()`, `estimateTagwiseDisp()` (`edgeR`) - フィルタリング: `filterByExpr()` \| - 発現変動解析: `exactTest()`, `glmFit()`, `glmLRT()` (`edgeR`) - 正規化: `calcNormFactors()` (`edgeR`) \|
MSnSet（Mass Spectrometry Set）	Bioconductor	質量分析データ（プロテオミクスデータ）の統合管理。	`MSnbase`, `MSstats`	- `exprs()`, `pData()`, `fData()` を用いた抽出	- データ操作: `exprs()`, `pData()`, `fData()` (`MSnbase`) - データ整形: `normalize()`, `logTransform()`	- 発現変動解析: `msqrob()`, `diffExp()` (`MSstats`) - 正規化: `normalize()` (`MSnbase`) \|
xcmsSet	Bioconductor	質量分析データのピーク検出とアライメント。	`xcms`, `CAMERA`	- `peaks()`, `group()` を用いた抽出	- ピーク検出: `findPeaks.centWave()`, `findPeaks.massifquant()` (`xcms`) - データ整形: `group()`, `fillPeaks()` (`xcms`)	-発現変動解析: `diffreport()` (`CAMERA`) - ピークアノテーション: `annotate()` (`CAMERA`)

データ構造の理解が高度な解析の基礎

各解析において取り扱うデータの複雑性に応じて、非常に多様なデータ構造が定義されている。
実際には、パッケージ特有のデータ構造が定義されている場合も多い。
あるパッケージの関数を基に解析フローを構築する場合、以下のスキルが必要になる:
- 各関数の入力として適切なデータ構造に合わせる。
- 自らデータ構造を構築する、または別の構造から変換する操作を行う。
出力結果を別のパッケージの関数に入力する際は、以下が必要:
- 出力データから必要な情報を抽出。
- データ構造を適切に調べ、操作する。
高度なデータ解析における基礎は、データ構造を的確に理解し、自在に操作できるスキルである。
データ構造の理解と操作が解析の基盤となることを強調し、Rにおけるデータ操作の本質を探る。

基本的なデータ構造と操作

Rにおける基本的なデータ構造は、ベクトル、行列、データフレーム、リストである。

ベクトル（Vector）

Rの最も基本的なデータ構造。
数値、文字列、論理値などを1次元で格納する。
c() 関数を用いて作成し、要素の追加、削除、並べ替えが可能。

例:

vec <- c(1, 2, 3, 4, 5)
vec[3]  # 3番目の要素を取得
vec[vec > 2]  # 2より大きい要素を抽出

行列（Matrix）

同じデータ型の要素を2次元で格納するデータ構造。
matrix() 関数を用いて作成し、行列演算や部分行列の取得が可能。

例:

mat <- matrix(1:9, nrow = 3)
mat[2, 3]  # 2行3列目の要素を取得

データフレーム（Data Frame）

異なるデータ型の列を持つことができる2次元データ構造。
data.frame() 関数で作成し、列名や行名でアクセス可能。
表形式のデータを扱うのに最適。

例:

df <- data.frame(
  name = c("Alice", "Bob", "Charlie"),
  age = c(25, 30, 35),
  score = c(80, 90, 85)
)
df$name  # "name" 列を取得

リスト（List）

異なるデータ型や構造を持つ要素を格納できるデータ構造。
list() 関数を用いて作成し、リスト内の要素には名前やインデックスでアクセスできる。

例:

lst <- list(name = "Alice", age = 25, scores = c(80, 90, 85))
lst$scores  # ベクトル "scores" を取得

データ構造を理解し自在に操作する

データ構造を理解するだけでは不十分であり、それらは相補的に使用される。
例えば、行列は数値データを効率的に格納するが、行列の操作にはベクトルが頻繁に使われる。
以下に、行列とベクトルの操作例を示す。

行列とベクトルの関係

行列とベクトルは、Rにおけるデータ操作の基本単位。
これらを組み合わせることで効率的なデータ処理が可能。

例1: 行列の作成とベクトルでの列抽出

# 行列の作成（2行3列の行列）
matrix_data <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)
# 行列の表示
print(matrix_data)

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

この行列は、2行3列で構成されている。

# 列を指定するベクトル
columns_to_extract <- c(2, 3)
# ベクトルを使って行列から特定の列を抽出
extracted_columns <- matrix_data[, columns_to_extract]
# 抽出した列の表示
print(extracted_columns)

     [,1] [,2]
[1,]    3    5
[2,]    4    6

ベクトルで列番号を指定して、行列から特定の列を抽出している。

例2: 行列の行ごとの和の計算

# 行ごとの和を計算
row_sums <- rowSums(matrix_data)
# 行ごとの和の表示
print(row_sums)

[1]  9 12

rowSums() 関数を使って、各行の要素を合計している。

例3: 行列の要素にベクトルを用いた条件付き処理

# 元の行列を表示
print(matrix_data)

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

# 条件を満たす要素に1を加算
matrix_data[matrix_data >= 3] <- matrix_data[matrix_data >= 3] + 1
# 結果を表示
print(matrix_data)

     [,1] [,2] [,3]
[1,]    1    4    6
[2,]    2    5    7

行列の条件付き操作を行い、特定の要素に1を加算している。

行列の構造を操作するということは、ベクトルで操作するということ

行列の操作は、内部的には「ベクトル」の操作である。

行列とベクトルの関係

Rでは、行列は単一のベクトルとして管理されており、列方向に要素が並べられたベクトル。
行列操作はベクトル操作の延長として理解できる。

# 行列の作成（2行3列の行列）
matrix_data <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)
print(matrix_data)

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

この行列は内部的にベクトル c(1, 2, 3, 4, 5, 6) として格納されている。

条件付き操作におけるベクトルの役割

行列の条件付き操作もベクトルとして処理される。

# 条件に基づく論理ベクトルの生成
logical_vector <- matrix_data >= 3
print(logical_vector)

     [,1]  [,2]  [,3]
[1,] FALSE  TRUE  TRUE
[2,] FALSE  TRUE  TRUE

# 条件を満たす要素に1を加算
matrix_data[logical_vector] <- matrix_data[logical_vector] + 1
print(matrix_data)

     [,1] [,2] [,3]
[1,]    1    4    6
[2,]    2    5    7

まとめ

行列はベクトルの集合体:
- 行列は内部的にはベクトルとして保存されている。
- 行列の操作はベクトルの操作と同じように考えることができる。
条件付き処理もベクトルとして処理:
- 行列の条件付き操作もベクトルとして処理され、論理ベクトルが生成される。
ベクトル操作の延長としての行列操作:
- 行列操作はベクトル操作の延長として理解でき、柔軟で効率的なデータ操作が可能になる。

他の構造に対する操作もベクトルによって操作する

他のデータ構造の操作においても、行列構造と同様に、ベクトルによる操作が基本となる。
データフレームでは、列ごとに異なるデータ型を格納でき、ベクトルを用いて効率的な操作が可能である。

データフレームにおけるベクトル操作の応用

データフレームの各列は内部的にベクトルとして格納され、効率的にデータ操作が可能。

列の選択やフィルタリング

例1: 列の抽出

# データフレームの作成
df <- data.frame(
  ID = 1:5,
  Name = c("Alice", "Bob", "Charlie", "David", "Eve"),
  Score = c(85, 92, 78, 90, 87)
)
print(df)

  ID    Name Score
1  1   Alice    85
2  2     Bob    92
3  3 Charlie    78
4  4   David    90
5  5     Eve    87

# Name列とScore列を抽出
name_vector <- df[, c(2, 3)]
print(name_vector)

      Name Score
1    Alice    85
2      Bob    92
3  Charlie    78
4    David    90
5      Eve    87

df データフレームの Name と Score 列を抽出し、ベクトル name_vector として格納。

例2: 条件に基づく行のフィルタリング

# Scoreが85以上の行を抽出
filtered_df <- df[df$Score >= 85, ]
print(filtered_df)

  ID   Name Score
1  1  Alice    85
2  2    Bob    92
4  4  David    90
5  5    Eve    87

Score 列が85以上の行を条件付きでフィルタリング。

リストにおけるベクトル操作の応用

リストは、ベクトルを拡張した構造で、異なるデータ型や構造をまとめて保持できる柔軟なデータ構造。

例1: ベクトルによる要素の選択

# 数値ベクトルを格納したリストの作成
num_list <- list(a = 1:5, b = 6:10, c = 11:15)
print(num_list)

$a
[1] 1 2 3 4 5

$b
[1]  6  7  8  9 10

$c
[1] 11 12 13 14 15

# 選択する要素の名前をベクトルとして定義
select_vector <- c("a", "c")

# リストから選択した要素を抽出
selected_elements <- num_list[select_vector]
print(selected_elements)

$a
[1] 1 2 3 4 5

$c
[1] 11 12 13 14 15

リスト num_list から a と c の要素をベクトル select_vector を用いて抽出。

例2: リスト内の各要素に対して操作を行う

# リスト内の各ベクトルに対して 2 を加算
modified_list <- lapply(num_list, function(x) x + 2)
print(modified_list)

$a
[1] 3 4 5 6 7

$b
[1]  8  9 10 11 12

$c
[1] 13 14 15 16 17

lapply 関数を使って、リスト内の各ベクトルに対して +2 の操作を行う。

例3: 条件付きフィルタリング

# 各ベクトルから、条件を満たす要素のみを抽出
filtered_list <- lapply(num_list, function(x) x[x > 8])
print(filtered_list)

$a
integer(0)

$b
[1]  9 10

$c
[1] 11 12 13 14 15

lapply を使い、リスト内の各ベクトルに対して、条件 x > 8 を満たす要素のみを抽出。

シングルセルデータにおけるベクトル操作の応用

SingleCellExperiment オブジェクトは、シングルセルRNA-seqデータを管理するためのデータ構造であり、細胞や遺伝子に関する情報を統合的に管理する。
ベクトル操作は、このオブジェクトに対するフィルタリングやデータ抽出の際に役立つ。

シングルセルデータの条件付きフィルタリング

例1: 細胞数のフィルタリング

library(SingleCellExperiment)

# SingleCellExperiment オブジェクトの作成（ダミーデータ）
sce <- SingleCellExperiment(
  assays = list(counts = matrix(rpois(100, lambda = 5), nrow = 10, ncol = 10)),
  colData = DataFrame(CellType = rep(c("A", "B"), each = 5))
)

# 各細胞の総カウント数を計算
total_counts <- colSums(counts(sce))

# 総カウント数が 20 以上の細胞を選択
sce_filtered <- sce[, total_counts >= 20]
print(dim(sce_filtered))

[1] 10  9

各細胞の総カウント数を計算し、20以上の細胞をフィルタリングする。
total_counts >= 20 の条件ベクトルを使って細胞をサブセット化。

ベクトル操作のまとめ

Rのデータ構造におけるベクトル操作は、基本的かつ重要な役割を果たす。
ベクトルはデータのサブセット化、フィルタリング、演算などに使用され、複雑な解析の基盤となる。
主要なデータ構造に対して普遍的に使えるベクトル操作を以下にまとめる。

1. ベクトル基本操作

ベクトルの生成: c(1, 2, 3, 4, 5)、1:5、seq(1, 5, by = 1)
要素のアクセス: x[1]（1番目の要素）
条件付き抽出: x[x > 3]（3より大きい要素）
要素の追加: c(x, 6)（末尾に要素を追加）
演算: x + 1（各要素に1を加算）

2. データフレームにおけるベクトル操作

列の抽出: df$column または df[["column"]]
行の抽出: df[1, ]（1行目）
条件付きフィルタリング: df[df$column > 10, ]
列の演算: df$new_column <- df$column1 + df$column2

3. 行列におけるベクトル操作

行列の要素: matrix_data[1, ]（1行目）
条件付き置換: matrix_data[matrix_data > 5] <- 0
行列積: matrix_data %*% c(1, 2, 3)
列の抽出: matrix_data[, 2]（2列目）

4. リストにおけるベクトル操作

リストの要素抽出: list_data[[1]]
ベクトルでフィルタリング: lapply(list_data, function(x) x[x > 3])
リストの要素長取得: sapply(list_data, length)
リスト内のベクトル演算: lapply(list_data, function(x) x + 1)

5. SummarizedExperimentやSingleCellExperimentにおけるベクトル操作

遺伝子やサンプルのフィルタリング: se[rowData(se)$gene == "TP53", ]
特定サンプルの抽出: se[, colData(se)$condition == "treated"]
発現データの操作: assay(se)["gene1", ] <- assay(se)["gene1", ] + 1
アノテーション情報の操作: rowData(se)$type <- "oncogene"

6. DGEListにおけるベクトル操作

条件ごとの抽出: dge$samples[group == "control", ]
発現データのフィルタリング: dge$counts[dge$counts > 5]
サンプルの正規化: dge$samples$norm.factors <- calcNormFactors(dge)

シークエンサーの定量データを管理するためのオブジェクト構造

BioconductorのSummarizedExperimentやSingleCellExperimentは、バイオインフォマティクス解析において非常に有用なデータ構造である。
定量値、メタ情報、アノテーション情報を一つのオブジェクト内で統一的に管理し、解析の整合性と効率を向上させる。
フィルタリングやサブセット化を行う際に、関連情報が自動的に更新されることでデータの一貫性が保たれる。

主要な利点

データの整合性を保つ
- Assayデータ、メタデータ、アノテーションデータを一元管理し、データの整合性を確保する。
柔軟なデータ操作
- 必要に応じてアノテーションやサンプル情報を追加・調整できる。
他の解析パッケージとの連携
- DESeq2やedgeRなどの解析パッケージと連携が可能であり、解析パイプラインを効率化。

補足： 前身のExperimentSet - ExperimentSetはかつてのデータ構造で、GEOqueryパッケージではまだ使用されることがあるが、SummarizedExperimentに変換することで解析の柔軟性を高める。

SummarizedExperimentオブジェクトの概要

図1: SummarizedExperimentオブジェクトの構造

Assays: 遺伝子発現量などの実データが含まれる行列（行が遺伝子、列がサンプル）。
Row Data: 各遺伝子に対応するアノテーション情報（例: 遺伝子名、染色体情報）。
Column Data: 各サンプルに対応するメタデータ（例: サンプルの状態、年齢）。

SummarizedExperimentオブジェクトの作成

library(SummarizedExperiment)

# 発現データ
assay_data <- matrix(rnorm(1000), nrow=100, ncol=10)

# サンプルメタデータ
col_data <- DataFrame(
  condition = rep(c("Control", "Treatment"), each=5),
  age = c(30, 32, 29, 45, 50, 28, 36, 33, 44, 40)
)

# 遺伝子アノテーションデータ
row_data <- DataFrame(
  gene_name = paste0("Gene", 1:100),
  chromosome = rep(c("chr1", "chr2", "chr3"), length.out=100)
)

# SummarizedExperimentオブジェクトの作成
se <- SummarizedExperiment(assays = list(counts = assay_data), 
                           rowData = row_data, 
                           colData = col_data)

# 結果の確認
se

SummarizedExperimentの基本操作

Assayデータの操作

Assayデータの取得

# Assayデータを取得
assay_data <- assay(se)
head(assay_data)

Assayデータの変更

# Assayデータの変更
assay(se)[1, ] <- rnorm(10)

Row Dataの操作

Row Dataの取得

# Row Dataを取得
row_data <- rowData(se)
head(row_data)

Row Dataの変更

# 遺伝子アノテーションの変更
rowData(se)$gene_name[1] <- "NewGene1"

Column Dataの操作

Column Dataの取得

# Column Dataを取得
col_data <- colData(se)
head(col_data)

Column Dataの変更

# サンプルメタデータの変更
colData(se)$condition[1] <- "NewCondition"

エラー処理とデータ構造の確認

データ構造を確認し、問題がないかをチェックする。

# データ構造を確認
str(se)

データが適切に読み込まれているか、データ型が正しいかを確認し、必要に応じて修正する。

`str()` 関数の結果を読み解く

str() 関数は、データ構造の概要を簡潔に表示するための強力なツールであり、特に複雑なデータ構造を持つ SummarizedExperiment や SingleCellExperiment のようなオブジェクトの内部構造を理解する際に役立つ。

`str()` 関数で表示される主な情報：

データのクラス: オブジェクトがどのクラスに属しているかを表示（例：SummarizedExperiment, SingleCellExperiment）。
オブジェクトの長さと構造: 各要素の数や階層構造（データフレームであれば行数と列数、リストであれば各要素の個数）。
各要素の型と内容: 各要素のデータ型（数値、文字列、因子、リストなど）と内容。

`SummarizedExperiment` オブジェクトの `str()` 結果の読み解き方

SummarizedExperiment は、発現データ（アッセイデータ）、遺伝子情報（rowData）、サンプル情報（colData）を統合管理する。

出力例

Formal class 'SummarizedExperiment' [package "SummarizedExperiment"] with 5 slots
  ..@ assays    :Formal class 'Assays' [package "SummarizedExperiment"] with 1 slot
  .. ..@ data:List of 1
  .. .. ..$ counts: int [1:20, 1:5] 15 8 7 14 12 10 11 7 16 12 ...
  ..@ rowData   :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  ..@ colData   :Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
  .. .. ..$ condition: Factor w/ 3 levels "A","B","C": 1 1 2 2 3
  ..@ metadata  :list()
  ..@ NAMES     :NULL
  ..@ elementMetadata: NULL

出力例の解釈

Formal class ‘SummarizedExperiment’: このオブジェクトは、SummarizedExperiment クラスに属している。
@ assays: 発現データを格納するスロット。list 型のデータに、counts 行列が含まれており、20行5列の整数値（int 型）が入っている。
@ rowData: 遺伝子情報（rowData）を格納するスロット。今回は詳細は表示されていないが、通常は遺伝子アノテーションが格納される。
@ colData: サンプルに関するメタデータを格納するスロット。ここでは condition という因子（Factor）変数があり、3つのレベル（“A”, “B”, “C”）を持っている。
@ metadata: メタデータを格納するリストで、今回は空である。
@ NAMES, elementMetadata: 特定の名前や要素のメタデータが格納されるスロットで、今回は NULL である。

`SummarizedExperiment`データ変更の方法について

SummarizedExperiment のデータ変更は、assays、rowData、colData といったスロットを操作して行う。

1. Assayデータの変更

Assayデータの取得と変更

# Assayデータの取得
assay_data <- assay(se)

# 例: 最初の遺伝子の発現データを全てゼロに変更
assay(se)[1, ] <- 0

新たなAssayデータの追加

# 正規化された新たなアッセイデータを追加
normalized_counts <- log2(assay_data + 1)
assays(se)$normalized_counts <- normalized_counts

2. Row Dataの変更

Row Dataの取得と変更

# Row Dataの取得
row_data <- rowData(se)

# 例: 最初の遺伝子名を "NewGene1" に変更
rowData(se)$gene_name[1] <- "NewGene1"

新たなRow Dataの追加

# 新しいアノテーションデータ（例: パスウェイ情報）を追加
rowData(se)$pathway <- c("Pathway1", "Pathway2", "Pathway3", ..., "Pathway100")

3. Column Dataの変更

Column Dataの取得と変更

# Column Dataの取得
col_data <- colData(se)

# 例: 最初のサンプルの条件を "NewCondition" に変更
colData(se)$condition[1] <- "NewCondition"

新たなColumn Dataの追加

# 新しいメタデータ（例: サンプルの処理時間）を追加
colData(se)$treatment_time <- c(24, 48, 72, 24, 48)

4. データの削除

特定の遺伝子やサンプルを削除することで、不要な情報を取り除ける。

# 例: 最初の遺伝子を削除
se <- se[-1, ]

# 例: 最初のサンプルを削除
se <- se[, -1]

`SingleCellExperiment` オブジェクトの構造と活用

SingleCellExperiment オブジェクトは、SummarizedExperiment の基盤の上に作られており、シングルセルRNAシーケンスデータを管理・解析するために設計されたデータ構造。
以下の3つの主要なコンポーネントを持つ：

Assays: 遺伝子発現量などのデータが含まれる行列。行が遺伝子や特徴量、列が細胞に対応。
Row Data: 各遺伝子に関連するアノテーション情報（遺伝子名や機能、アノテーションなど）。
Column Data: 各細胞に対応するメタデータ（細胞の状態やクラスター情報など）。

次元削減結果は、reducedDims スロットに追加されることがある（PCA、UMAPなど）。

SingleCellExperiment オブジェクトの作成

シングルセルRNAシーケンスデータを基に SingleCellExperiment オブジェクトを作成する方法。

library(SingleCellExperiment)

# 発現データ (行が遺伝子、列が細胞)
assay_data <- matrix(rnorm(1000), nrow=100, ncol=10)

# 細胞に関するメタデータ
col_data <- DataFrame(
  cluster = rep(c("Cluster1", "Cluster2"), each=5),
  condition = rep(c("Control", "Treatment"), each=5)
)

# 遺伝子に関するアノテーションデータ
row_data <- DataFrame(
  gene_name = paste0("Gene", 1:100),
  chromosome = rep(c("chr1", "chr2", "chr3"), length.out=100)
)

# SingleCellExperiment オブジェクトを作成
sce <- SingleCellExperiment(assays = list(counts = assay_data), 
                            rowData = row_data, 
                            colData = col_data)

基本的な操作

SummarizedExperiment に類似した操作が可能で、次元削減データの操作が追加されている。

Assay データの操作

Assay データの取得

assay_data <- assay(sce)
head(assay_data)

Assay データの変更

assay(sce)[1, ] <- rnorm(10)

Row Data の操作

Row Data の取得

row_data <- rowData(sce)
head(row_data)

Row Data の変更

rowData(sce)$gene_name[1] <- "NewGene1"

Column Data の操作

Column Data の取得

col_data <- colData(sce)
head(col_data)

Column Data の変更

colData(sce)$condition[1] <- "NewCondition"

次元削減データの操作

次元削減データの追加

pca_result <- prcomp(t(assay(sce)))
reducedDims(sce)$PCA <- pca_result$x

次元削減データの取得

umap_result <- reducedDims(sce)$UMAP

str() 関数を用いたデータ構造の確認

SingleCellExperiment オブジェクトの内部構造を確認するために str() を使用。

str(sce)

出力例の解釈：

Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 6 slots
  ..@ int_elementMetadata: Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
  ..@ colData            : Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
  ..@ reducedDims        : Formal class 'SimpleList' [package "S4Vectors"] with 1 slot
  ..@ assays             :Formal class 'Assays' [package "SummarizedExperiment"] with 1 slot
  .. ..@ data:List of 1
  .. .. ..$ counts: int [1:10, 1:10] 14 10 15 14 5 12 16 16 5 7 ...
  ..@ metadata           : list()
  ..@ int_colData        : Formal class 'DataFrame' [package "S4Vectors"] with 6 slots

データの変更方法

1. Assayデータの変更

Assayデータの取得と変更

assay_data <- assay(sce)
assay(sce)[1, ] <- 0

新たなAssayデータの追加

normalized_counts <- log2(assay_data + 1)
assays(sce)$normalized_counts <- normalized_counts

2. Row Dataの変更

Row Dataの取得と変更

rowData(sce)$gene_name[1] <- "NewGene1"

新たなRow Dataの追加

rowData(sce)$pathway <- c("Pathway1", "Pathway2", "Pathway3", ..., "Pathway100")

3. Column Dataの変更

Column Dataの取得と変更

colData(sce)$condition[1] <- "NewCondition"

新たなColumn Dataの追加

colData(sce)$treatment_time <- c(24, 48, 72, 24, 48)

4. Reduced Dimensionの操作

次元削減データの追加

reducedDims(sce)$PCA <- pca_result$x

次元削減データの取得

umap_result <- reducedDims(sce)$UMAP

5. データの削除

特定の行（遺伝子）や列（サンプル）の削除：

sce <- sce[-1, ]  # 最初の遺伝子を削除
sce <- sce[, -1]  # 最初のサンプルを削除

6. データの確認とエラー処理

データ構造の確認

str(sce)

データの結合

Rにおけるデータ操作で重要な概念の一つはデータの結合。これは異なるデータセットを一貫性を保ちながら結合し、分析を進めるために必要な操作。
データの結合には、行方向の結合、列方向の結合、およびキーによる結合といった方法がある。

1. 行方向の結合（縦結合）

行方向の結合は、同じ列構造を持つ二つ以上のデータフレームや行列を、行単位で連結する操作。
異なる実験条件やサンプルのデータを一つのデータセットにまとめる際に利用される。

関数の例

rbind(): 同じ列名を持つデータフレームを結合。
- 例: rbind(df1, df2)
bind_rows() (dplyr パッケージ): 異なる列名を持つデータも柔軟に結合。
- 例: bind_rows(df1, df2)

使用例

# データフレームの作成
df1 <- data.frame(ID = 1:3, Name = c("Alice", "Bob", "Charlie"))
df2 <- data.frame(ID = 4:6, Name = c("David", "Eve", "Frank"))
# 行方向の結合
df_combined <- rbind(df1, df2)
print(df_combined)

2. 列方向の結合（横結合）

列方向の結合は、同じ行数を持つデータフレームや行列を列単位で結合する操作。
同じサンプルに異なる変数を追加する際に使用。

関数の例

cbind(): 行数が一致するデータを列方向に結合。
- 例: cbind(df1, df2)
bind_cols() (dplyr パッケージ): 列名の衝突を自動で解決。
- 例: bind_cols(df1, df2)

使用例

# データフレームの作成
df1 <- data.frame(ID = 1:3, Score = c(85, 90, 88))
df2 <- data.frame(Name = c("Alice", "Bob", "Charlie"), Age = c(25, 30, 28))
# 列方向の結合
df_combined <- cbind(df1, df2)
print(df_combined)

3. キーによる結合（データベース結合）

キーによる結合は、異なるデータセット間の共通のキーを基にデータを結びつける操作。SQLのJOINと類似。
異なるデータソースから情報を統合して新しいデータセットを作成。

関数の例

merge(): 共通の列（キー）を基にデータフレームを結合。
- 例: merge(df1, df2, by = "ID")
left_join(), inner_join(), full_join() (dplyr パッケージ): 左結合、内部結合、全結合などを柔軟に実施。
- 例: left_join(df1, df2, by = "ID")

使用例

# データフレームの作成
df1 <- data.frame(ID = 1:3, Name = c("Alice", "Bob", "Charlie"))
df2 <- data.frame(ID = 2:4, Score = c(92, 85, 78))
# キーによる結合
df_combined <- merge(df1, df2, by = "ID")
print(df_combined)

結合タイプごとの説明

1. 左結合（Left Join）

左側のデータセットを基準に、右側のデータセットを結合。
左側に存在するが右側に存在しない場合はNA。

left_join(df1, df2, by = "ID")

  ID    Name Score Subject Grade
1  1   Alice    85    <NA>  <NA>
2  2     Bob    92    Math     A
3  3 Charlie    78 Science     B
4  4   David    90    <NA>  <NA>

2. 内部結合（Inner Join）

両方のデータセットに共通のキーが存在する行のみを結合。

inner_join(df1, df2, by = "ID")

  ID    Name Score Subject Grade
1  2     Bob    92    Math     A
2  3 Charlie    78 Science     B

3. 全結合（Full Join）

両方のデータセットすべての行を保持し、欠損部分は NA。

full_join(df1, df2, by = "ID")

  ID    Name Score Subject Grade
1  1   Alice    85    <NA>  <NA>
2  2     Bob    92    Math     A
3  3 Charlie    78 Science     B
4  4   David    90    <NA>  <NA>
5  5    <NA>  <NA> History     C

4. 右結合（Right Join）

右側のデータセットを基準に結合し、左側に一致しない行はNA。

right_join(df1, df2, by = "ID")

  ID    Name Score
1  2     Bob    92
2  3 Charlie    85
3  4    <NA>    78

5. アンチ結合（Anti Join）

左側に存在し、右側に存在しない行のみを抽出。

anti_join(df1, df2, by = "ID")

  ID    Name Score
1  1   Alice    85
2  4   David    90

4. 結合における注意点

共通キーの有無: データセット間で共通のキーがない場合、結合が正しく行えない。
欠損値の扱い: 外部結合では欠損値が発生するため、事前に考慮が必要。
列名の重複: 列名が重複する場合、結合前に調整するか確認が必要。

演習：SummarizedExperimentとSingleCellExperimentにおけるデータ結合操作

SummarizedExperiment と SingleCellExperiment を使ったデータ結合操作を学ぶ。
外部アノテーションデータの統合や、複数のシングルセルデータの結合を実施。

演習1: HGNCデータベースからのアノテーション追加（SummarizedExperiment）

SummarizedExperiment オブジェクトを作成し、仮のHGNCデータを使用して遺伝子シンボルを追加。
left_join() 関数を使って、SummarizedExperiment の rowData に遺伝子シンボルを統合。

library(SummarizedExperiment)

# サンプルデータを作成
counts <- matrix(rpois(100, lambda = 10), ncol = 5)
rownames(counts) <- paste0("gene", 1:20)  # 仮の遺伝子ID
colData <- DataFrame(condition = c("A", "A", "B", "B", "C"))
se <- SummarizedExperiment(assays = list(counts = counts), colData = colData,)

# 外部アノテーションデータ
hgnc_data <- data.frame(
  GeneID = paste0("gene", 1:20),  
  GeneSymbol = paste0("Symbol", 1:20),
  stringsAsFactors = FALSE
)

# rowDataにアノテーションを追加
row_data <- as.data.frame(rowData(se))  
joined_data <- left_join(row_data, hgnc_data, by = c("rowname" = "GeneID"))
rowData(se) <- DataFrame(joined_data)

# 結果を表示
print(rowData(se))

確認事項:

left_join() でオリジナルの遺伝子情報を保持しつつ、アノテーションが正しく追加されたか確認。
他の結合方法（inner_join, full_join, right_join, merge）の比較を行う。

演習2: Ensemblデータベースからのアノテーション追加（SummarizedExperiment）

biomaRt パッケージを使用してEnsemblから遺伝子シンボル情報を取得し、SummarizedExperiment オブジェクトに統合。
getBM() 関数と left_join() によるアノテーション統合プロセスを学ぶ。

# 必要なライブラリをロード
library(SummarizedExperiment)
library(dplyr)
library(biomaRt)

# サンプルデータ作成 (実際のENSG IDを使うため、ランダムなものではなく現実のものを使用)
counts <- matrix(rpois(100, lambda = 10), ncol = 5)
# 例として最初の20個の現実のENSG IDを用いる
ensembl_gene_ids <- c(
  "ENSG00000139618", "ENSG00000241860", "ENSG00000227232", "ENSG00000157764",
  "ENSG00000142192", "ENSG00000243989", "ENSG00000198786", "ENSG00000240361",
  "ENSG00000163435", "ENSG00000198947", "ENSG00000160712", "ENSG00000112715",
  "ENSG00000141510", "ENSG00000171862", "ENSG00000135486", "ENSG00000119866",
  "ENSG00000137757", "ENSG00000166913", "ENSG00000165699", "ENSG00000134460"
)
rownames(counts) <- ensembl_gene_ids
colData <- DataFrame(condition = c("A", "A", "B", "B", "C"))
se <- SummarizedExperiment(assays = list(counts = counts), colData = colData)

# Ensemblデータベースからアノテーションを取得
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
annotation_data <- getBM(
  attributes = c("ensembl_gene_id", "hgnc_symbol"),
  filters = "ensembl_gene_id",
  values = rownames(counts),
  mart = ensembl
)

# アノテーションをSummarizedExperimentに追加
row_data_df <- as.data.frame(rowData(se))
# 遺伝子情報の結合を行う
merged_data <- merge(data.frame(rownames = rownames(se)), annotation_data, by.x = "rownames", by.y = "ensembl_gene_id", all.x = TRUE)
rowData(se) <- DataFrame(merged_data)

# 結果を表示
print(rowData(se))

確認事項:

遺伝子シンボル情報が正しく統合されているか確認。

演習3: シングルセルデータのバッチ結合（SingleCellExperiment）

2つの異なるバッチのシングルセルデータを作成し、cbind() を使ってデータを結合。
SingleCellExperiment の結合が正しく行われるか確認。

library(SingleCellExperiment)

# シングルセルRNA-seqデータの作成
expr_data1 <- matrix(rpois(20, lambda = 5), nrow = 5)
expr_data2 <- matrix(rpois(20, lambda = 5), nrow = 5)

# バッチ情報を含むサンプルアノテーションを作成
batch_info1 <- DataFrame(SampleID = paste0("Cell", 1:4), Batch = "Batch1")
batch_info2 <- DataFrame(SampleID = paste0("Cell", 5:8), Batch = "Batch2")

# SingleCellExperimentオブジェクトの作成
sce1 <- SingleCellExperiment(assays = list(counts = expr_data1), colData = batch_info1)
sce2 <- SingleCellExperiment(assays = list(counts = expr_data2), colData = batch_info2)

# バッチ間のデータ結合
combined_sce <- cbind(sce1, sce2)

# サンプル情報と発現データを確認
print(colData(combined_sce))
print(assay(combined_sce, "counts"))

確認事項:

バッチ間でシングルセルデータが正しく結合され、colDataにバッチ情報が保持されているか確認。

データの変換・要約

概要: データ変換・要約処理は、特徴の把握、フィルタリング、可視化に重要。総和やZスコア変換、分散・平均の計算などが行われる。
バイオインフォマティクスでの応用:
- フィルタリングのスコア計算
- 可視化用の統計データ作成

フィルタリングの役割と重要性

目的: 重要でない特徴を除外し、データ軽量化や偽陽性の排除を行う。
例:
- 遺伝子発現の総和、平均、分散・標準偏差を基に不要なデータを除外。
- 効率的な解析やモデル構築の支援。

可視化と統計データの利用

役割: データの探索的調査や可視化のカスタマイズが必要。
品質管理:
- 外れ値確認: Zスコア基づくヒートマップ、低次元空間の散布図。
- バッチ効果確認: 散布図、ボックスプロット、密度分布による可視化。
使用ツール: ggplot2, pheatmap などを利用しつつ、自らのデータ変換とカスタム可視化が求められる。

効率的なデータ操作と要約の重要性

一般的な処理: データフレームや行列に対する処理・要約。
理解するべき点: 表形式のデータ構造に対する操作。

使用する主要な関数とパッケージ

apply() 関数と dplyr パッケージ:
- apply() は行列やデータフレームに対して関数を行や列ごとに適用。
- dplyr は効率的なデータ操作をサポート。

1. apply() 関数

用途: 行列やデータフレームに関数を適用するために使用。
例1: 行ごとの合計計算

# 行列データ作成
data <- matrix(rnorm(20), nrow = 5)

# 行ごとの合計を計算
row_sums <- apply(data, 1, sum)
print(row_sums)

例2: 列ごとの平均計算

# 列ごとの平均を計算
col_means <- apply(data, 2, mean)
print(col_means)

1 は行方向、2 は列方向を示す。関数をデータ構造に適用。

applyファミリーの他の関数

lapply(): リスト要素ごとに関数を適用。
sapply(): リストやベクトルに関数を適用し、結果をベクトルや行列で返す。

関数名	説明	使い方の例	主な用途と有用な場面
apply()	行列やデータフレームに対して、行や列ごとに関数を適用する。	`apply(data, 1, sum)` で行ごとの合計を計算	行列の各行または各列に同じ処理を施す場合。例: 行列の行ごとの合計、列ごとの平均を計算する。
lapply()	リストやベクトルに対して関数を適用し、リストとして結果を返す。	`lapply(my_list, mean)` でリスト内の平均を計算	リストの各要素に対して同じ処理を行い、その結果をリスト形式で保持したい場合。例: リスト内の各ベクトルに対して集計処理を行う。
sapply()	`lapply()` に似ているが、結果を可能であればベクトルや行列として返す。	`sapply(my_list, mean)` でベクトルとして結果を返す	`lapply()` と同様の操作を行いつつ、より簡潔な形式（ベクトルや行列）で結果を取得したい場合。例: データフレームの各列に対して平均を計算し、ベクトルとして返す。
vapply()	`sapply()` に似ているが、返される値の型を指定できる。	`vapply(my_list, mean, numeric(1))` で型を指定	`sapply()` のように結果をベクトルや行列として返しつつ、型を保証したい場合。例: 処理の結果が特定の型（数値、文字列）であることを確実にしたい場合。
mapply()	複数のベクトルやリストに対して、要素ごとに関数を適用する。	`mapply(sum, 1:5, 6:10)` で要素ごとの合計を計算	複数の引数を持つ関数に対して、複数のベクトルやリストを並行して処理したい場合。例: 2つのベクトルの要素ごとに処理を行う。
tapply()	ベクトルの値をグループごとに分割し、各グループに関数を適用する。	`tapply(scores, group, mean)` でグループごとの平均を計算	データをグループ化して、そのグループごとに関数を適用したい場合。例: サンプルをグループごとに分け、その中で平均や合計を計算する。
rapply()	リストのネストされた要素に対して再帰的に関数を適用する。	`rapply(nested_list, sqrt)` でネスト内の要素に平方根を適用	ネストされたリストに対して再帰的に処理を行いたい場合。例: 複雑なリスト構造内のすべての数値要素に対して操作を行う。
eapply()	環境に含まれる全てのオブジェクトに関数を適用する。	`eapply(env, mean)` で環境内のオブジェクトに処理を適用	環境内のすべてのオブジェクトに対して関数を適用したい場合。例: 環境内のデータフレームに対して一括処理を行う。

applyファミリーで実装可能な処理
- 行・列の総和、最大値、z-scoreなどは apply() を使って実装可能
- 行列やデータフレームに対して効率的な集計や変換を行うために、専用の関数が提供されている
- 特に大規模なデータ処理では、各遺伝子やサンプルごとの特徴量計算やフィルタリングが重要

行方向・列方向の操作に使える関数一覧

関数	説明	使用例
`rowSums()`	各行の合計値を計算する	`rowSums(matrix_data)`
`colSums()`	各列の合計値を計算する	`colSums(matrix_data)`
`rowMeans()`	各行の平均値を計算する	`rowMeans(matrix_data)`
`colMeans()`	各列の平均値を計算する	`colMeans(matrix_data)`
`scale()`	行列やデータフレームの列を標準化（Zスコア）する	`scale(matrix_data)`
`t()`	行列の転置を行う（行と列を入れ替える）	`t(matrix_data)`
`rowMin()`	各行の最小値を計算する（`matrixStats` パッケージが必要）	`rowMin(matrix_data)`
`colMin()`	各列の最小値を計算する（`matrixStats` パッケージが必要）	`colMin(matrix_data)`
`rowMax()`	各行の最大値を計算する（`matrixStats` パッケージが必要）	`rowMax(matrix_data)`
`colMax()`	各列の最大値を計算する（`matrixStats` パッケージが必要）	`colMax(matrix_data)`
`rowVars()`	各行の分散を計算する（`matrixStats` パッケージが必要）	`rowVars(matrix_data)`
`colVars()`	各列の分散を計算する（`matrixStats` パッケージが必要）	`colVars(matrix_data)`
`rowMedians()`	各行の中央値を計算する（`matrixStats` パッケージが必要）	`rowMedians(matrix_data)`
`colMedians()`	各列の中央値を計算する（`matrixStats` パッケージが必要）	`colMedians(matrix_data)`

使用場面の例

rowSums() / colSums():
- 遺伝子発現データに対して、各遺伝子ごと（行）の発現総量や各サンプルごと（列）の発現総量を計算する際に便利。
scale():
- データを標準化して解析に適したスケールに変換する際に使用（Zスコア変換など）。
matrixStats 関連の関数:
- matrixStats パッケージには、行列に対する分散や中央値、最小値、最大値などの計算を効率的に行う関数が豊富に用意されているため、特に大規模な行列データに対して使用することが多い。

使用例

行列データの標準化
- 行列データを標準化する際、Zスコア変換を行うことがよくある。
- apply() 関数や scale() 関数を使って実装できる。

例3: 列ごとにZスコアを計算

# 行列データを標準化
scaled_data <- scale(data)
print(scaled_data)

scale() 関数は、各列の平均を0、標準偏差を1に標準化する。
- 異なるスケールのデータを統一されたスケールに変換可能。

行ごと・列ごとの特徴量の抽出とフィルタリング
- 変換や要約の目的でフィルタリングを行うことが重要。
- 分散が低い遺伝子をフィルタリングして除外したり、特定の条件を満たすサンプルを抽出する。

例4: 分散が低い遺伝子をフィルタリング

# 各行の分散を計算
variances <- apply(data, 1, var)

# 分散がしきい値以上の行を選択
filtered_data <- data[variances > 0.5, ]
print(filtered_data)

apply() 関数で各行の分散を計算。
- 分散が一定以上の遺伝子だけを抽出。

複数の変換ステップの適用
- フィルタリング後に標準化するなど、複数の変換ステップを組み合わせることが一般的。

例5: フィルタリングと標準化の組み合わせ

# 分散フィルタリング後、標準化
filtered_scaled_data <- scale(filtered_data)
print(filtered_scaled_data)

フィルタリングされたデータを標準化して、解析に適したスケールに変換する。

演習：`apply()` ファミリーおよび行・列方向の操作に使える関数

演習1: 行ごとの合計を計算する

次の行列データ data に対して、apply() 関数を使用して行ごとの合計を計算してみよう。

data <- matrix(1:9, nrow = 3)

# 行ごとの合計を計算
row_sums <- apply(data, 1, sum)
print(row_sums)

演習2: 列ごとの平均を計算する

問題
上記の data 行列に対して、列ごとの平均を計算するコードを実装しなさい。apply() または colMeans() 関数を使用すること。

# apply() を使用した列ごとの平均
col_means <- apply(data, 2, mean)
print(col_means)

# colMeans() を使用した列ごとの平均
col_means_alt <- colMeans(data)
print(col_means_alt)

演習3: 標準偏差の計算

上記の行列に対して、各列の標準偏差を計算するコードを実装してみよう。apply() または matrixStats パッケージの colSds() を使用する。

# apply() を使用した標準偏差の計算
col_sds <- apply(data, 2, sd)
print(col_sds)

# matrixStats パッケージを使用した標準偏差の計算
library(matrixStats)
col_sds_alt <- colSds(data)
print(col_sds_alt)

演習4: 分散の計算とフィルタリング

次のデータセット gene_expression において、行（遺伝子ごと）の分散を計算し、分散が 0.2 以上の遺伝子だけをフィルタリングしてみよう。

gene_expression <- matrix(runif(50, 0, 1), nrow = 10)

# 各行の分散を計算
variances <- apply(gene_expression, 1, var)

# 分散がしきい値以上の行をフィルタリング
filtered_data <- gene_expression[variances > 0.2, ]
print(filtered_data)

演習5: 行列の標準化（Zスコア変換）

gene_expression データに対して、列ごとのZスコアを計算し、各遺伝子（行）ごとの標準化されたスコアを計算してみよう。scale() 関数を使用する。

# 列ごとの標準化（Zスコア変換）
scaled_data <- scale(gene_expression)
print(scaled_data)

演習6: apply() ファミリーの活用

次のリスト my_list に対して、lapply() 関数を使用して、各ベクトルの合計を計算してみよう。

my_list <- list(a = 1:5, b = 6:10, c = 11:15)

# 各ベクトルの合計を計算
list_sums <- lapply(my_list, sum)
print(list_sums)

演習7: 列ごとの最大値を計算する

上記の gene_expression データに対して、列ごとの最大値を計算してみよう。apply() または matrixStats パッケージの colMaxs() 関数を使用する。

# apply() を使用した列ごとの最大値の計算
col_maxs <- apply(gene_expression, 2, max)
print(col_maxs)

# matrixStats パッケージを使用した列ごとの最大値の計算
library(matrixStats)
col_maxs_alt <- colMaxs(gene_expression)
print(col_maxs_alt)

演習8: 条件付きフィルタリング

次のデータセット data_filter において、各行の平均が0.5以上の行だけをフィルタリングしてみよう。

data_filter <- matrix(runif(50, 0, 1), nrow = 10)

# 各行の平均を計算
row_means <- apply(data_filter, 1, mean)

# 平均が0.5以上の行をフィルタリング
filtered_data <- data_filter[row_means >= 0.5, ]
print(filtered_data)

演習9: 列方向の最小値の計算

次のデータセット data_stats において、各列の最小値を計算してみよう。apply() 関数または matrixStats パッケージの colMins() 関数を使う。

data_stats <- matrix(rnorm(100), nrow = 10)

# apply() を使用した列ごとの最小値の計算
col_mins <- apply(data_stats, 2, min)
print(col_mins)

# matrixStats パッケージを使用した列ごとの最小値の計算
library(matrixStats)
col_mins_alt <- colMins(data_stats)
print(col_mins_alt)

2. dplyr パッケージによる要約操作

dplyr パッケージでは、summarise() や mutate() を使って複雑な変換や要約処理が可能。
group_by()を組み合わせることで、条件ごとの集約や変換も効率的に処理できる。
データを可読性の高い形式で要約し、複雑な条件付き処理にも対応可能。

例6: 列ごとの特徴量を計算

library(dplyr)

# サンプルデータフレームを作成
df <- data.frame(
  ID = 1:5,
  Score = c(85, 92, 78, 90, 87),
  Group = c("A", "B", "A", "B", "A")
)

# グループごとの平均スコアを計算
mean_score <- df %>%
  group_by(Group) %>%
  summarise(mean = mean(Score))
print(mean_score)

group_by() と summarise() を使い、Score 列を Group でグループ化し、それぞれのグループの平均を計算。
複数の条件付き処理に対応し、効率的に計算できる点が優れている。

例7: 新しい列を追加してデータを変換

# Zスコアを計算して新しい列を追加
df <- df %>%
  mutate(Z_Score = (Score - mean(Score)) / sd(Score))
print(df)

mutate() を使って Score のZスコアを計算し、新たな列としてデータフレームに追加。
変換や新しい特徴量の追加が簡潔かつ可読性の高いコードで実現。

例8: 条件付きのデータ変換と要約

# Group "A" のみにZスコアを計算して追加
df <- df %>%
  group_by(Group) %>%
  mutate(Z_Score = if_else(Group == "A", (Score - mean(Score)) / sd(Score), NA_real_))
print(df)

mutate() と if_else() を組み合わせて、Group が “A” の場合のみ Zスコアを計算。
条件付きの変換処理もシンプルに記述でき、複雑な処理において優位性がある。

3. long/wide型のデータ形式変換

ggplot2などの可視化パッケージでデータをプロットする際や、複数カテゴリのデータを比較する場合に、表データの形式変換（long型⇔wide型）が必要となる場面が多い。
dplyr と tidyr パッケージを組み合わせることで、データ形式変換が容易になる。

long型とwide型のデータ形式変換の例

wide型データ: 遺伝子発現量や疾患状態、年齢、性別といった情報を含むデータを例とする。

wide型データ例:

SampleID	GeneA	GeneB	GeneC	Disease_Status	Age	Sex
1	5.4	3.2	7.1	Healthy	45	Male
2	6.1	2.8	6.9	Disease	60	Female
3	5.9	3.5	7.4	Disease	50	Male
4	5.2	3.1	7.0	Healthy	38	Female

wide型からlong型への変換
- pivot_longer()を使用して遺伝子発現情報を一列にまとめ、Gene列に遺伝子名を、Expression列に発現量を格納。

library(dplyr)
library(tidyr)

# サンプルデータ
df <- data.frame(
  SampleID = c(1, 2, 3, 4),
  GeneA = c(5.4, 6.1, 5.9, 5.2),
  GeneB = c(3.2, 2.8, 3.5, 3.1),
  GeneC = c(7.1, 6.9, 7.4, 7.0),
  Disease_Status = c("Healthy", "Disease", "Disease", "Healthy"),
  Age = c(45, 60, 50, 38),
  Sex = c("Male", "Female", "Male", "Female")
)

# wide型からlong型への変換
long_df <- df %>%
  pivot_longer(cols = starts_with("Gene"), 
               names_to = "Gene", 
               values_to = "Expression")

print(long_df)

結果（long型データ）:

SampleID	Disease_Status	Age	Sex	Gene	Expression
1	Healthy	45	Male	GeneA	5.4
1	Healthy	45	Male	GeneB	3.2
1	Healthy	45	Male	GeneC	7.1
2	Disease	60	Female	GeneA	6.1
2	Disease	60	Female	GeneB	2.8
2	Disease	60	Female	GeneC	6.9

long型からwide型への変換
- pivot_wider()を使ってGene列を展開し、wide型に戻す。

# long型からwide型へ変換
wide_df <- long_df %>%
  pivot_wider(names_from = Gene, values_from = Expression)

print(wide_df)

結果（wide型データ）:

SampleID	Disease_Status	Age	Sex	GeneA	GeneB	GeneC
1	Healthy	45	Male	5.4	3.2	7.1
2	Disease	60	Female	6.1	2.8	6.9
3	Disease	50	Male	5.9	3.5	7.4
4	Healthy	38	Female	5.2	3.1	7.0

long/wide形式変換の意義

long型: 各サンプルの複数の遺伝子発現を1列にまとめることで、ggplot2などでカテゴリ別に遺伝子発現量をプロットする際に便利。
wide型: データの確認や要約統計を計算する際に適している。

applyファミリーとdplyrの比較

applyファミリー: 行列やリストに対する単純な集計や変換に強みがある。
dplyr: 複数の変換ステップや条件付き処理に優れ、コードの可読性が高い。データフレーム操作やフィルタリング、グループ化に強力。

使い分けの基準

行列やリストの操作 → applyファミリー
データフレーム操作、グループ化や条件付き処理 → dplyr

複雑なオブジェクトに対する処理について

SummarizedExperiment や SingleCellExperiment オブジェクトでは、データ構造から発現データやアノテーション情報に対するデータ変換や要約処理が可能。
これらのオブジェクト自体は、発現データとアノテーションの管理を行うためのデータ構造であり、ユーザーが基本関数や他のパッケージの関数を利用して変換・要約処理を行う。
特定の解析における変換・要約処理は、DESeq2 や edgeR のようなBioconductorのパッケージで自動的に処理される。

発現データの行列に対する変換処理

発現データ（assays スロットに格納されている行列）は通常の行列操作が可能。
遺伝子ごとの分散や平均、Zスコアへの変換など、基本的な統計的処理を行うことができる。

例: 行ごとの平均と分散を計算

# SummarizedExperimentオブジェクトの発現データを取得
counts <- assay(se)

# 遺伝子ごとの平均と分散を計算
gene_means <- rowMeans(counts)
gene_vars <- apply(counts, 1, var)

print(gene_means)
print(gene_vars)

assay() 関数を使って発現データを取得し、rowMeans() で遺伝子ごとの平均、apply() で遺伝子ごとの分散を計算。
通常の行列操作を SummarizedExperiment の発現データに対してそのまま適用可能。

可視化

可視化は、データ解析において重要な技術であり、単なる形式や見栄えだけではなく、データを正確に伝え、理解を深めるためのプロセスである。
データを視覚化する際には、見る人が正確なメッセージを受け取れるようにデザインし、意図的な偏りや情報過多による誤解を避けることが大切である。
良い可視化の基準は、視覚的に明確なメッセージを伝え、適切な情報を強調することである。可視化の効果を最大限に引き出すためには、目的に応じた適切な可視化手法とツールを選ぶ必要がある。

1. 可視化の目的と役割

可視化の目的の多様性

可視化には次の3つの主な目的がある。

データの傾向を直感的に理解するため（探索的解析）
- 初期段階のデータ解析で、パターンや異常値を発見し、仮説を立てるために行う。
- 例: ggplot2 や base plot を使用したシンプルで迅速な可視化。
仮説を検証し、重要な情報を強調するため（説明的解析）
- グループ間の比較や相関関係を示す際に、特定のデータを強調する。
- 例: 棒グラフ、散布図、箱ひげ図を使用して、仮説を視覚的に裏付ける。
発表・公表用に高品質な可視化を作成するため
- 論文やプレゼン用の可視化では、美しさや分かりやすさが重要。
- 例: ggplot2 や plotly を使用して、視覚的に美しい図を作成。

下流解析における具体的役割

技術的バイアスやバッチ効果によるグループ構造の検出
- 可視化を用いて、技術的要因による影響を検出する。
- 例: PCAやUMAPプロットでバッチごとの分離を確認。
生物学的なグループ構造の発見
- データ内の生物学的な特徴やグループ構造を明らかにする。
- 例: UMAPプロットで疾患状態ごとのサンプルグループを発見。
品質管理
- 異常値やノイズを検出し、データが適切かどうかを確認する。
- 例: ボックスプロットやヒートマップで異常なサンプルや遺伝子を検出。
各処理・解析ステップごとの方法の評価
- フィルタリングや正規化などの解析ステップの結果を可視化して評価する。
- 例: 正規化前後の分布をヒストグラムで比較。
仮説の確からしさ（実験・介入効果の可視化）
- 実験条件による発現変動を視覚化し、仮説を確認する。
- 例: ボルケーノプロットで有意な発現変動遺伝子を確認。
分子間の相互作用やモジュールの発見
- 分子間のネットワークを可視化し、重要な分子やモジュールを発見する。
- 例: ネットワーク解析で分子相互作用を可視化。
データ全体の俯瞰
- 大量データを一望できる可視化を行い、解析の方向性を定める。
- 例: ヒートマップで発現パターンの似た遺伝子群を俯瞰。

演習1: 技術的バイアスやバッチ効果の検出

目的: PCAプロットを用いて、技術的バイアス（バッチ効果）がサンプル間にあるかどうかを確認する。

# パッケージをロード
library(ggplot2)
library(DESeq2)

# サンプルのデータ作成 (例: RNA-seqデータ)
set.seed(123)
counts <- matrix(rpois(10000, lambda = 10), ncol = 10)
batch <- factor(rep(c("Batch1", "Batch2"), each = 5))
colData <- data.frame(batch = batch)

# DESeq2オブジェクトを作成
dds <- DESeqDataSetFromMatrix(countData = counts, colData = colData, design = ~batch)

# 正規化
vsd <- vst(dds)

# PCAプロット
pcaData <- plotPCA(vsd, intgroup = "batch", returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color = batch)) +
  geom_point(size = 3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  ggtitle("PCA plot with batch effect")

確認事項: サンプルがバッチごとに分かれてクラスタリングされているかどうかを確認し、技術的バイアスが存在するかを検討する。

演習2: 生物学的なグループ構造の発見

目的: UMAPプロットを用いて、生物学的なグループ（例: 健康 vs. 疾患）を発見する。

# パッケージをロード
library(ggplot2)
library(umap)

# サンプルデータ作成
set.seed(123)
data <- matrix(rnorm(2000), nrow = 200, ncol = 10)
group <- factor(rep(c("Control", "Disease"), each = 5))

# UMAPを実行
umap_result <- umap(data)

# 結果をプロット
umap_df <- data.frame(UMAP1 = umap_result$layout[,1],
                      UMAP2 = umap_result$layout[,2],
                      group = group)

ggplot(umap_df, aes(UMAP1, UMAP2, color = group)) +
  geom_point(size = 3) +
  ggtitle("UMAP plot showing biological groups")

確認事項: UMAPプロットでサンプルが健康（Control）と疾患（Disease）のグループで明確に分かれているかを確認する。

演習3: 品質管理のための可視化

目的: ボックスプロット、PCAに基づく主成分スコアの散布図、ヒートマップを用いて、サンプル間で異常なデータポイントがないかを確認する。

# パッケージをロード
library(ggplot2)
library(DESeq2)

# サンプルデータ作成（通常サンプルと異常サンプルを含む）
set.seed(123)
counts <- matrix(rpois(900, lambda = 100), ncol = 9)  # 通常サンプル
abnormal_sample <- rpois(100, lambda = 500)  # 異常サンプル（外れ値）
counts <- cbind(counts, abnormal_sample)  # 異常サンプルを追加
colnames(counts) <- paste0("Sample", 1:10)
colData <- data.frame(condition = factor(rep(c("A", "B"), each = 5)))

# DESeq2オブジェクトを作成
dds <- DESeqDataSetFromMatrix(countData = counts, colData = colData, design = ~ condition)

# ボックスプロットで各サンプルのリード数の分布、バイアスを確認
boxplot(counts, main = "Boxplot with Abnormal Sample for Quality Control", 
        xlab = "Samples", ylab = "Counts", outline = TRUE)

# varianceStabilizingTransformation() を使用して正規化
vsd <- varianceStabilizingTransformation(dds)

# PCAプロットを作成
pcaData <- plotPCA(vsd, intgroup = "condition", returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color = condition)) +
  geom_point(size = 3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  ggtitle("PCA Plot with Abnormal Sample for Quality Control")

確認事項: 異常に高いまたは低いサンプルがないか、サンプル間でバイアスがないかを確認する。

演習4: 各解析ステップの評価

目的: 正規化前後のRNA-seqデータのヒストグラムを比較し、正規化の効果を確認する。

# パッケージをロード
library(ggplot2)
library(DESeq2)

# サンプルデータ作成 (RNA-seqデータ)
set.seed(123)
counts <- matrix(rpois(1000, lambda = 10), ncol = 10)

# DESeq2オブジェクトを作成
dds <- DESeqDataSetFromMatrix(countData = counts, colData = data.frame(condition = factor(rep(c("A", "B"), each = 5))), design = ~ condition)

# 正規化前のヒストグラム
par(mfrow = c(1, 2))
hist(log2(counts + 1), breaks = 50, main = "Before normalization", xlab = "Log2(counts + 1)")

# varianceStabilizingTransformation()を使用して正規化
vsd <- varianceStabilizingTransformation(dds)

# 正規化後のヒストグラム
hist(assay(vsd), breaks = 50, main = "After normalization", xlab = "VST counts")

確認事項: 正規化前後の分布を確認し、正規化によりデータが整ったかどうかを判断する。

演習5: 仮説の確からしさ（発現変動の可視化）

目的: ボルケーノプロットで発現変動遺伝子を視覚化し、実験の結果を確認する。

# パッケージをロード
library(ggplot2)
library(DESeq2)

# サンプルデータ作成
set.seed(123)
dds <- makeExampleDESeqDataSet()

# 発現変動解析
dds <- DESeq(dds)
res <- results(dds)

# 欠損値の処理
res$log2FoldChange[is.na(res$log2FoldChange)] <- 0
res$pvalue[is.na(res$pvalue)] <- 1

# 遺伝子名をラベルとして追加
res$gene <- rownames(res)

# ボルケーノプロット
ggplot(res, aes(x = log2FoldChange, y = -log10(pvalue))) +
  geom_point(alpha = 0.5) +
  ggtitle("Volcano Plot") +
  xlab("Log2 Fold Change") +
  ylab("-Log10 P-value") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "red") +
  
  # 閾値を満たす遺伝子にラベルを追加
  geom_text(
    data = subset(res, abs(log2FoldChange) > 2 & pvalue < 0.05), 
    aes(label = gene), 
    vjust = 1.5, 
    color = "blue", 
    size = 3
  )

確認事項: 有意に発現変動している遺伝子がボルケーノプロット上で強調されているか確認する。

演習6: 分子間の相互作用の可視化

目的: ネットワーク解析に基づく分子間の相互作用を可視化する。

# パッケージをロード
library(igraph)

# サンプルネットワークデータ作成
nodes <- data.frame(name = paste0("Gene", 1:10))
edges <- data.frame(from = sample(nodes$name, 10, replace = TRUE),
                    to = sample(nodes$name, 10, replace = TRUE))

# グラフオブジェクトの作成
graph <- graph_from_data_frame(edges, vertices = nodes, directed = FALSE)

# ネットワークプロット
plot(graph, vertex.label = V(graph)$name, main = "Gene Interaction Network")

確認事項: 分子間の相互作用を視覚化し、重要なノードや相互作用を発見できるか確認する。

演習7: データ全体の俯瞰

目的: ヒートマップを使用して、発現パターンを俯瞰し、遺伝子やサンプルの全体像を把握する。

# パッケージをロード
library(pheatmap)

# サンプルデータ作成
set.seed(123)
data <- matrix(rnorm(100), nrow = 10)

# ヒートマップを作成
pheatmap(data, cluster_rows = TRUE, cluster_cols = TRUE, 
         main = "Heatmap of Gene Expression")

2. 可視化の基本的な要素

データの構造に基づく可視化の選択

適切なグラフを選択することで、視覚的に情報を正確に伝える。データの構造や伝えたい情報に応じて異なるグラフを使用する。
1. 1変数データ（単一変数データ）
  
  データの分布や要約統計量を視覚的に示す。
  - ヒストグラム
    
    データの分布を直感的に理解するためのグラフ。
    - 使用例: 学生のテストスコアの分布。
  - 箱ひげ図（Box Plot）
    
    中央値、四分位範囲、最小値・最大値、外れ値を視覚化。
    - 使用例: 遺伝子発現データの分布の比較。
2. 2変数データ（2つの変数の関係）
  
  2つの変数の関係性を視覚化。
  - 散布図（Scatter Plot）
    
    変数間の相関やトレンドを視覚化。
    - 使用例: 遺伝子発現量と治療効果の相関。
  - バイオリンプロット
    
    分布を詳細に示すためのプロット。
    - 使用例: 条件間の遺伝子発現の比較。
  - 箱ひげ図（Box Plot）
    
    外れ値や中央値の比較が容易。
    - 使用例: サンプル間でのRNA-seqデータの分布の比較。
3. 多変数データ（高次元データ）
  
  複数の変数を同時に扱うための可視化手法。
  - ペアプロット（Pair Plot）
    
    すべての変数の組み合わせを可視化。
    - 使用例: 複数の遺伝子発現データの相関関係。
  - ファセットによる可視化
    
    データをカテゴリごとに分割して比較。
    - 使用例: 複数サンプルにおける遺伝子発現の分布。
  - 次元削減後の可視化（PCA, UMAP, t-SNE）
    
    高次元データを2次元や3次元に圧縮。
    - 使用例: シングルセルRNA-seqデータのクラスタリング。
ラベルや軸の設計

軸ラベルやタイトルは、データの意味を明確に伝えるために重要。スケールや目盛りも適切に設定する必要がある。
- ラベルの明確さ
  
  軸ラベルやタイトルは簡潔で情報を明確にするべき。
  - 例: 「遺伝子発現量 (Log2)」「サンプルグループ」
- 軸のスケールと目盛り
  
  データの範囲に応じたスケールを選択し、目盛りを読みやすく配置。
カラースキームの選択

色はデータの情報を強調するために使用。
- 見やすさと情報伝達のバランス
  
  カテゴリ変数や連続変数の表現に適した色を選ぶ。
  - 例: コントラストのある色を選択（例: 青と赤）。
- カラーブラインド対応
  
  カラーブラインド対応の配色（例: viridis）を使用。

3. コードによる実装例: 探索的解析と説明的解析

迅速なデータ探索のためのシンプルな可視化

plot(), hist(), boxplot() などを使用したベースRでの可視化。
高度な可視化パッケージによる詳細なカスタマイズ

ggplot2を使用して、軸ラベル、色分け、凡例、ファセットなどをカスタマイズ可能なグラフ作成。

迅速なデータ探索のためのシンプルな可視化
ベースRによる可視化の例
- plot(): 変数間の関係を確認するための散布図など、基本的なグラフを作成可能。
- hist(): データの分布を確認するヒストグラムを作成。
- boxplot(): サンプル間の分布を比較し、外れ値を特定するボックスプロットを描画。

例1: ベースRを使ったヒストグラムと散布図の作成

# サンプルデータ作成
set.seed(123)
data <- rnorm(100)

# ヒストグラム
hist(data, main = "Histogram of Data", xlab = "Value", col = "lightblue")

# 散布図
x <- rnorm(100)
y <- 2*x + rnorm(100)
plot(x, y, main = "Scatter Plot", xlab = "X", ylab = "Y")

ヒストグラム: データの分布が視覚化され、中心傾向や広がりが把握できる。

グラフ名	対象となる構造	使用場面	関数名	入力データの構造・形式	主な引数（データ、スケール、ラベル、軸、色、タイトル）
散布図	2変数（連続変数）	変数間の関係性を確認する	`plot()`	ベクトル、データフレーム（wide）	`x, y`（データ）, `xlim, ylim`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）
ヒストグラム	1変数（連続変数）	データの分布を確認する	`hist()`	ベクトル	`x`（データ）, `breaks`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）
ボックスプロット	1変数または2変数（カテゴリ）	データの範囲、中央値、外れ値を確認する	`boxplot()`	ベクトル、データフレーム（long）	`x, y`（データ）, `ylim`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）
バープロット	カテゴリデータ	カテゴリごとのカウントを視覚化	`barplot()`	テーブル、行列	`height`（データ）, `ylim`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）
ペアプロット	多変数（連続変数）	変数間の関係性を一度に確認	`pairs()`	データフレーム（wide）	`x`（データ）, `xlim, ylim`（スケール）, `labels`（ラベル）, `pch`（色）, `main`（タイトル）
線グラフ	2変数（連続変数、時間変数）	時系列データのトレンドを確認	`plot()`	ベクトル	`x, y`（データ）, `xlim, ylim`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）
パイチャート	カテゴリデータ	カテゴリごとの割合を視覚化	`pie()`	ベクトル、テーブル	`x`（データ）, `labels`（ラベル）, `col`（色）, `main`（タイトル）
ヒートマップ	行列（多変数）	行列形式のデータのパターンを視覚化	`heatmap()`	行列	`x`（データ）, `scale`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）
棒グラフ	カテゴリデータ	カテゴリごとの数値比較	`barplot()`	数値ベクトル、行列	`height`（データ）, `ylim`（スケール）, `xlab, ylab`（軸）, `col`（色）, `main`（タイトル）

散布図: 2つの変数の関係が視覚化され、相関の有無や線形関係が確認できる。

高度な可視化パッケージによる詳細なカスタマイズ

ggplot2を中心に説明。ggplot2は、高度にカスタマイズ可能な可視化ツールを提供し、発表や論文に使える品質のグラフを作成可能。
レイヤーベースのアプローチに基づき、軸ラベルや色分け、凡例の位置などを細かく調整できる。

1. ggplotにおける可視化の考え方

ggplot2 は Grammar of Graphics に基づいてデータを可視化。
文法（Grammar）は、可視化を構成する基本要素を組み合わせて、さまざまなグラフを作成可能。

Grammar of Graphicsの基本要素

データ（Data）
- グラフに表示する元となるデータセット。
- ggplot()の第一引数に指定。
審美的マッピング（Aesthetic Mapping, aes）
- データ変数を視覚要素（位置、色、大きさなど）に対応させる。
- 例: aes(x = variable1, y = variable2, color = group)
幾何オブジェクト（Geometries）
- グラフの形状を定義。例えば、散布図や棒グラフなど。
- 例: geom_point()（散布図）、geom_bar()（棒グラフ）
スケール（Scales）
- 審美的マッピングのデータをどのようにスケールするかを設定（軸や色の範囲など）。
- 例: scale_x_continuous(limits = c(0, 10))（x軸のスケール）
統計変換（Statistics）
- データの集計や変換を指定。例: トレンド線や要約統計量。
- 例: geom_smooth(method = "lm")（線形回帰線を追加）
座標系（Coordinate Systems）
- 描画に使用する座標系。通常は直交座標系だが、極座標なども可能。
- 例: coord_polar()（極座標で描画）
ファセット（Facets）
- グループごとのデータを分割して描画。
- 例: facet_wrap(~ group)（グループごとにサブプロットを作成）

`ggplot2`での可視化の流れ

基本構文
- ggplot2では、複数の関数を+で繋げて可視化を行う。

ggplot(data, aes(x = ..., y = ...)) +
  geom_<plot_type>(aes(...)) +
  labs(title = "Plot Title", x = "X-axis label", y = "Y-axis label") +
  theme(...)

可視化プロセス
1. データ形式の変換
  - long形式またはwide形式にデータを整形する。
  - pivot_longer() や pivot_wider() を使用。
```
dataset <- pivot_longer(wide_dataset)
dataset <- pivot_wider(long_dataset)
```
1. データを指定
  - 描画するデータセットを指定する。
```
ggplot(data = dataset)
```
1. 美学マッピングを指定
  - x軸、y軸、色などに対応する変数を指定する。
```
ggplot(data = dataset, aes(x = var1, y = var2, color = group))
```
1. 幾何オブジェクトを指定
  - 散布図、棒グラフなどの形式を指定。
```
+ geom_point()
```
1. スケールやテーマを設定
  - 軸範囲やカラースケール、テーマを調整する。
```
+ scale_x_continuous(limits = c(0, 10))
+ theme_minimal()
```
1. タイトルやラベルを設定
  - タイトルや軸ラベルを追加する。
```
+ ggtitle("Title")
+ xlab("X-axis Label")
+ ylab("Y-axis Label")
```

例1: 散布図を作成する流れ

# ggplotによる散布図の作成例
library(ggplot2)

# サンプルデータ作成
data <- data.frame(x = rnorm(100), y = rnorm(100), group = rep(1:2, each = 50))

# 散布図を作成
ggplot(data, aes(x = x, y = y, color = as.factor(group))) +  # データと美学マッピング
  geom_point() +                                             # 幾何オブジェクト
  ggtitle("Scatter Plot Example") +                          # タイトル
  xlab("X-axis") + ylab("Y-axis") +                         # ラベル
  theme_minimal()                                            # スケール・テーマ

ポイント
- xとyの変数を散布図の軸に対応させ、groupで色分け。
- geom_point()を使用して点をプロットし、タイトルと軸ラベルを追加。

long形式の重要性

long形式のデータは、ggplot2で異なる美学マッピング（aesthetic mapping）を行う際に重要。
例: 異なる条件（時間、治療、実験条件）ごとに数値を比較する場合や、複数の変数に対して色や形を割り当てる場合。
wide形式は、計算や統計解析に適しているが、ggplot2での可視化には適さない場合がある。
解析結果の例
- 発現変動解析などでは、p値や対数比変動、統計的有意性の情報が整理された形式で提供され、long形式に変換しやすい。
メタ情報・アノテーション情報
- サンプルや実験条件に関する情報（メタ情報）や、遺伝子シンボルやGOカテゴリなどのアノテーション情報も、可視化で重要な役割を果たす。
適切なデータ変換
- 解析オブジェクトや結果からデータを取り出し、適切にフィルタリングしてlong形式へ変換することにより、ggplot2での効果的な可視化が可能。

`ggplot2`による詳細な可視化一覧表

グラフ名	対象となる構造	使用場面	関数名	入力データの構造・形式	主な引数（データ、スケール、ラベル、軸、色、タイトル）
散布図	2変数（連続変数）	変数間の関係性を視覚的に表現	`geom_point()`	データフレーム（long）	`aes(x, y)`（データ）, `scale_x_continuous(), scale_y_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(color)`（色）, `ggtitle()`（タイトル）
ヒストグラム	1変数（連続変数）	データの分布を視覚的に表現	`geom_histogram()`	データフレーム（long）	`aes(x)`（データ）, `scale_x_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(fill)`（色）, `ggtitle()`（タイトル）
ボックスプロット	1変数または2変数（カテゴリ）	カテゴリ間の分布の比較を行う	`geom_boxplot()`	データフレーム（long）	`aes(x, y)`（データ）, `scale_y_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(fill)`（色）, `ggtitle()`（タイトル）
バープロット	カテゴリデータ	カテゴリごとの比較	`geom_bar()`	データフレーム（long）	`aes(x)`（データ）, `scale_y_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(fill)`（色）, `ggtitle()`（タイトル）
ペアプロット	多変数（連続変数）	変数間の関係性をペアで確認	`GGally::ggpairs()`	データフレーム（wide）	`aes()`（データ）, `axis.ticks`（スケール）, `labs()`（軸）, `aes(color)`（色）, `ggtitle()`（タイトル）
線グラフ	2変数（連続変数、時間変数）	時系列データや連続変数のトレンドを可視化	`geom_line()`	データフレーム（long）	`aes(x, y)`（データ）, `scale_x_continuous(), scale_y_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(color)`（色）, `ggtitle()`（タイトル）
バイオリンプロット	1変数または2変数（カテゴリ）	分布の形状を詳細に確認	`geom_violin()`	データフレーム（long）	`aes(x, y)`（データ）, `scale_y_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(fill)`（色）, `ggtitle()`（タイトル）
ヒートマップ	行列（多変数）	行列形式のデータを可視化	`geom_tile()`	データフレーム（long）	`aes(x, y, fill)`（データ）, `scale_fill_gradient()`（スケール）, `labs(x, y)`（軸）, `ggtitle()`（タイトル）
ファセットプロット	多変数	グループごとのプロット分割	`facet_wrap()`	データフレーム（long）	`aes(x, y)`（データ）, `facet_wrap(~variable)`（スケール）, `labs(x, y)`（軸）, `ggtitle()`（タイトル）
PCA・UMAPプロット	高次元データ	次元削減後のデータを可視化	`geom_point()`	データフレーム（wide/long）	`aes(x, y)`（データ）, `scale_x_continuous(), scale_y_continuous()`（スケール）, `labs(x, y)`（軸）, `aes(color)`（色）, `ggtitle()`（タイトル）

演習：`ggplot2`を使った実際の解析場面を想定した可視化

以下の演習では、SummarizedExperimentオブジェクトを基に、特定の遺伝子の発現量を可視化したり、発現変動解析の結果を可視化することを通じて、ggplot2を使ったさまざまな可視化手法を学ぶ。

演習1: 特定遺伝子の発現量を条件間で比較

目的: サンプルごとの特定遺伝子の発現量を、条件（例えば治療群とコントロール群）で比較するプロットを作成する。
手順: 特定遺伝子のデータを抽出し、ggplot2を使って箱ひげ図（ボックスプロット）を作成する。

# パッケージをロード
library(SummarizedExperiment)
library(ggplot2)

# サンプルデータ作成
set.seed(123)
counts <- matrix(rpois(1000, lambda = 10), ncol = 10)
rownames(counts) <- paste0("gene", 1:100)  # 遺伝子名を設定
colnames(counts) <- paste0("Sample", 1:10)  # サンプル名を設定

# colData（サンプル情報）の設定
colData <- DataFrame(condition = factor(rep(c("Control", "Treatment"), each = 5)))

# SummarizedExperimentオブジェクトを作成
se <- SummarizedExperiment(assays = list(counts = counts), colData = colData)

# 特定の遺伝子（例: gene1）の発現量を抽出
gene1_expr <- assay(se)["gene1", ]  # 遺伝子名 'gene1' を指定
gene1_data <- data.frame(Sample = colnames(se), Expression = gene1_expr, Condition = colData(se)$condition)

# 箱ひげ図を作成
ggplot(gene1_data, aes(x = Condition, y = Expression, fill = Condition)) +
  geom_boxplot() +
  ggtitle("Gene1 Expression by Condition") +
  xlab("Condition") +
  ylab("Expression Level") +
  theme_minimal()

演習2: 薬剤濃度依存の発現変化プロット

目的: 各サンプルにおける薬剤濃度に依存した発現変化を可視化する。
手順: 濃度ごとにサンプルをグループ化し、発現量の変化を折れ線グラフでプロットする。

# パッケージをロード
library(SummarizedExperiment)
library(ggplot2)

# サンプルデータ作成
set.seed(123)
counts <- matrix(rpois(1000, lambda = 10), ncol = 10)
rownames(counts) <- paste0("gene", 1:100)  # 遺伝子名を設定
colnames(counts) <- paste0("Sample", 1:10)  # サンプル名を設定

# colData（サンプル情報）の設定
concentration <- c(0, 0.1, 1, 10, 50, 0, 0.1, 1, 10, 50)
colData <- DataFrame(concentration = concentration)

# SummarizedExperimentオブジェクトを作成
se <- SummarizedExperiment(assays = list(counts = counts), colData = colData)

# 特定遺伝子の発現量と濃度データを抽出
gene_expr <- assay(se)["gene2", ]  # 遺伝子名 'gene2' を指定
gene_data <- data.frame(Sample = colnames(se), Expression = gene_expr, Concentration = colData(se)$concentration)

# 薬剤濃度依存の発現変化プロット
ggplot(gene_data, aes(x = Concentration, y = Expression)) +
  geom_point() +
  geom_line() +
  ggtitle("Gene2 Expression Change by Drug Concentration") +
  xlab("Concentration (µM)") +
  ylab("Expression Level") +
  theme_minimal()

演習3: 異常サンプルを散布図でラベル付き表示

目的: 発現データの中で異常なサンプルを散布図で可視化し、異常サンプルにラベルを付けて表示する。
手順: 散布図を作成し、閾値を超えた異常なサンプルにラベルを付けて表示する。

# サンプルデータ作成
set.seed(123)
counts <- matrix(rnorm(1000, mean = 10, sd = 2), ncol = 10)

# サンプル名を設定
colnames(counts) <- paste0("Sample", 1:10)

# サンプル条件データを作成
colData <- DataFrame(condition = factor(rep(c("Control", "Treatment"), each = 5)))

# SummarizedExperimentオブジェクトを作成
se <- SummarizedExperiment(assays = list(counts = counts), colData = colData)

# PCAを実行するためにデータを標準化
assay_data <- t(assay(se))  # サンプルを行方向に
pca_result <- prcomp(assay_data, scale. = TRUE)

# PCAの結果をデータフレームにまとめる
pca_data <- data.frame(Sample = colnames(se), 
                       PC1 = pca_result$x[, 1], 
                       PC2 = pca_result$x[, 2], 
                       Condition = colData(se)$condition)

# 異常サンプルを判断する閾値（例: PC1が2以上）
threshold <- 2

# 散布図に異常サンプルをラベルとして追加
ggplot(pca_data, aes(x = PC1, y = PC2, color = Condition)) +
  geom_point() +
  geom_text(aes(label = ifelse(abs(PC1) > threshold, Sample, "")), vjust = -1) +
  ggtitle("PCA Plot with Outliers Labeled") +
  xlab("PC1") +
  ylab("PC2") +
  theme_minimal()

演習4: 発現変動解析の結果からボルケーノプロットを作成

目的: 発現変動解析の結果テーブルからボルケーノプロットを作成する。
手順: DESeq2の解析結果をもとに、ボルケーノプロットを作成し、統計的に有意な遺伝子を強調する。

# パッケージをロード
library(DESeq2)
library(ggplot2)

# サンプルデータ作成
set.seed(123)
dds <- makeExampleDESeqDataSet()

# 発現変動解析を実行
dds <- DESeq(dds)
res <- results(dds)

# NAの補完
res$log2FoldChange[is.na(res$log2FoldChange)] <- 0
res$pvalue[is.na(res$pvalue)] <- 1

# ボルケーノプロットの作成
ggplot(res, aes(x = log2FoldChange, y = -log10(pvalue))) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "red") +
  ggtitle("Volcano Plot of Differential Gene Expression") +
  xlab("Log2 Fold Change") +
  ylab("-Log10 P-value") +
  theme_minimal()

演習5: wide形式からlong形式への変換と可視化

目的:

データがwide形式で保存されている場合、それをlong形式に変換してから可視化する。特定の遺伝子発現量をサンプルごとに比較し、条件間での発現量の違いを可視化する。

手順:

pivot_longer()関数を使って、wide形式のデータをlong形式に変換する。
ggplot2でボックスプロットを作成して、条件ごとの遺伝子発現量を比較する。

# パッケージをロード
library(tidyr)
library(ggplot2)

# サンプルデータ作成 (wide形式)
counts <- data.frame(
  Gene = paste0("Gene", 1:5),
  Sample1 = rpois(5, 10),
  Sample2 = rpois(5, 12),
  Sample3 = rpois(5, 14),
  Sample4 = rpois(5, 16),
  Sample5 = rpois(5, 18),
  Condition = factor(c("Control", "Control", "Treatment", "Treatment", "Control"))
)

# データを確認 (wide形式)
print(counts)

# wide形式からlong形式へ変換
long_counts <- pivot_longer(counts, cols = starts_with("Sample"), 
                            names_to = "Sample", values_to = "Expression")

# long形式データの確認
print(long_counts)

# 条件ごとのボックスプロットを作成
ggplot(long_counts, aes(x = Condition, y = Expression, fill = Condition)) +
  geom_boxplot() +
  ggtitle("Gene Expression by Condition") +
  xlab("Condition") +
  ylab("Expression Level") +
  theme_minimal()

演習6: 複数遺伝子の発現変動解析 (wideからlongへの変換)

目的:

複数の遺伝子について、条件間での発現量を比較するために、wide形式のデータをlong形式に変換し、条件ごとの発現量を可視化する。

手順:

wide形式の遺伝子発現データを用意。
pivot_longer()でデータをlong形式に変換。
ggplot2で条件ごとの発現量をボックスプロットで表示。

# パッケージをロード
library(tidyr)
library(ggplot2)

# サンプルデータ作成 (wide形式: 複数遺伝子の発現量)
rna_counts <- data.frame(
  Gene = paste0("Gene", 1:5),
  Control_Sample1 = rnorm(5, mean = 100, sd = 15),
  Control_Sample2 = rnorm(5, mean = 102, sd = 14),
  Treatment_Sample1 = rnorm(5, mean = 120, sd = 18),
  Treatment_Sample2 = rnorm(5, mean = 122, sd = 17)
)

# wide形式からlong形式へ変換
long_rna_counts <- pivot_longer(rna_counts, cols = -Gene,
                                names_to = c("Condition", "Sample"),
                                names_sep = "_", values_to = "Expression")

# 発現量の可視化
ggplot(long_rna_counts, aes(x = Condition, y = Expression, fill = Condition)) +
  geom_boxplot() +
  facet_wrap(~ Gene, scales = "free_y") +
  ggtitle("Gene Expression Comparison by Condition") +
  xlab("Condition") +
  ylab("Expression Level") +
  theme_minimal()

演習7: 時系列データのlong形式変換と可視化

目的:

時系列データをlong形式に変換し、各遺伝子の発現量が時間とともにどのように変化するかを可視化する。

手順:

wide形式で保持されている時間ごとの発現データをpivot_longer()でlong形式に変換。
ggplot2で各遺伝子の時系列プロットを作成。

# パッケージをロード
library(tidyr)
library(ggplot2)

# サンプルデータ作成 (wide形式: 時系列データ)
time_series <- data.frame(
  Gene = paste0("Gene", 1:5),
  Time0 = rnorm(5, mean = 100, sd = 10),
  Time1 = rnorm(5, mean = 105, sd = 12),
  Time2 = rnorm(5, mean = 110, sd = 14),
  Time3 = rnorm(5, mean = 115, sd = 13)
)

# wide形式からlong形式へ変換
long_time_series <- pivot_longer(time_series, cols = -Gene, 
                                 names_to = "Time", values_to = "Expression")

# 時系列プロットを作成
ggplot(long_time_series, aes(x = Time, y = Expression, group = Gene, color = Gene)) +
  geom_line() +
  geom_point() +
  ggtitle("Gene Expression Over Time") +
  xlab("Time") +
  ylab("Expression Level") +
  theme_minimal()

演習8: 遺伝子間の相関関係の可視化 (wideからlongへの変換)

目的:

複数の遺伝子間の相関を計算し、その相関行列をlong形式に変換してヒートマップを作成する。

手順:

遺伝子発現データから相関行列を計算。
melt()を用いて相関行列をlong形式に変換。
ggplot2でヒートマップを作成。

# パッケージをロード
library(tidyr)
library(ggplot2)
library(reshape2)

# サンプルデータ作成 (wide形式: 複数遺伝子の発現量)
gene_data <- data.frame(
  Gene1 = rnorm(10, mean = 100, sd = 15),
  Gene2 = rnorm(10, mean = 110, sd = 20),
  Gene3 = rnorm(10, mean = 120, sd = 25),
  Gene4 = rnorm(10, mean = 130, sd = 30)
)

# 遺伝子間の相関行列を計算
cor_matrix <- cor(gene_data)

# 相関行列をlong形式に変換
long_cor_matrix <- melt(cor_matrix)

# ヒートマップで相関を可視化
ggplot(long_cor_matrix, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0) +
  ggtitle("Gene Correlation Heatmap") +
  xlab("Gene") +
  ylab("Gene") +
  theme_minimal()

バイオインフォマティクス解析におけるRの実践的基礎

Yusuke Matsui

2024-10-10

なぜRをあえて「講義」で学ぶのか？

R vs Python: 二元論に意味はあるのか？

AI時代のR学習: その意義とは？

Rを学ぶ上で「講義形式」が重要な理由

講義形式で学ぶことの意義: 体系的な知識と実践的スキル

実際の解析シナリオでのRの活用法を学ぶ

講義の目的

Rがバイオインフォマティクスで活用される場面

Rの役割と他ツールとの連携の重要性

例：RNA-seq解析におけるRの具体的な役割

RNA-seqパイプラインにおける各処理ステップと対応ツール

各ツールとスクリプトの役割

Rの具体的な役割

Rがバイオインフォマティクスで活躍できる場面と限界のある場面

大規模データ処理におけるRの限界

Rの強みを最大限に活かす場面

適材適所でツールを使い分ける

解析フローの設計

R、シェル、Pythonを組み合わせる際のベストプラクティス

Rの役割は下流解析において重要

Rを用いてできることの全貌

Rを用いた下流解析

Bioconductorパッケージを使うとできること

Bioconductorの特筆すべき利点

R + Bioconductorの強力な組み合わせ

Pythonとの比較

Rを用いた下流解析で必要な操作スキルとは？

状況に対応するためのスキル

Rでの解析における「親切な関数」と「不親切な関数」の存在

関数の出力データの分析・整理の必要性

必要なスキルを習得するためのステップ

例1: RNA-seq発現変動解析（limma）

各ステップの説明

定量アッセイに対する下流解析の一般的フローと必要なR操作

1. データの統合とアノテーション

2. データの前処理

3. 次元削減とクラスタリング

4. 機械学習と予測モデル

5. 発現変動解析と機能的アノテーション

6. ネットワーク解析

Rでデータを自在に操るための基礎

1. データの整形と前処理

2. データのフィルタリングとサブセット化

3. データのマージとアノテーション

4. 関数の出力結果の整形と解析

5. カスタムの可視化

6. パイプラインの構築と自動化

パッケージ関数の活用を前提としたR操作の本質的な基礎とは？

1. 入力・出力データのファイル形式およびデータ形式の多様性

2. Rにおけるデータ構造の理解と操作スキル

3. データの整形とサマリー作成

4. 可視化のためのデータ整形と統計的要約

データのライフサイクルを理解する

バイオインフォマティクス解析におけるデータの多様性

データ読み込みの基盤スキル

1. データ形式の理解

2. 引数やオプションの柔軟な使用

3. 試行錯誤とデータの確認

4. 想定外のデータフォーマットへの対応力

5. データがどうしても読み込めない場合の対応

6. シェルスクリプトによる事前処理

1. head でデータの確認

2. sed や awk を用いたデータ整形

3. grep で特定行を抽出

7. データ品質チェックとクリーニング

データ読み込みスキルの演習例題

例題1: CSVファイルの読み込みと整形

課題:

Rでのデータ生成と保存

演習:

解答例:

解説：

ステップ 1: CSVファイルの読み込み

ステップ 2: GeneID列を文字列型にする

ステップ 3: データの確認

例題2: タブ区切りファイルの読み込みと前処理

課題:

1. `head` でデータの確認

2. `sed` や `awk` を用いたデータ整形

3. `grep` で特定行を抽出

ステップ 2: `GeneID`列を文字列型にする

`str()` 関数の結果を読み解く

`str()` 関数で表示される主な情報：