一つにまとまったページ
Content from スプレッドシートによるデータ整理
最終更新日:2026-11-07 | ページの編集
所要時間: 60分
概要
質問
- 表形式データの整理方法は?
目的
- スプレッドシートとその長所・短所について学びましょう。
- データを効果的に活用するために、スプレッドシートではどのようにデータを整形すればよいでしょうか。
- スプレッドシートでよくあるエラーと、その修正方法について学びましょう。
- 整然データの原則に従ってデータを整理してください。
- カンマ区切り(CSV)やタブ区切り(TSV)などのテキストベースのスプレッドシート形式について学びましょう。
このエピソードは、Data Carpentriesの「生態学者のためのRによるデータ分析と可視化」レッスンに基づいています。
スプレッドシートソフト
質問
- データを適切に整理するためのスプレッドシート活用の基本原則とは何でしょうか?
目的
- コンピュータがデータセットを最大限に活用できるよう、データを整理する際のベストプラクティスについて説明してください。
要点
- 適切なデータ整理は、あらゆる研究プロジェクトの基盤となります。
適切なデータ整理は、研究プロジェクトの基盤となります。ほとんどの研究者は、スプレッドシートでデータを扱ったり、データ入力を行ったりしています。スプレッドシートソフトは、データテーブルの設計や、ごく基本的なデータ品質管理機能を利用する上で非常に有用なグラフィカルインターフェースです。 Broman and Woo (2018). も参照してください。
スプレッドシートの構成
スプレッドシートはデータ入力に適しています。そのため、スプレッドシートには膨大なデータが蓄積されています。研究者としての時間の多くは、この「データ加工」の段階に費やされることになります。決して楽しい作業ではありませんが、不可欠なものです。ここでは、データの整理方法や、より効率的なデータ加工のための実践的な手法について解説します。
このレッスンで学べないこと
- スプレッドシートで統計を行う方法
- スプレッドシートでグラフ作成を行う方法
- スプレッドシートソフトでコードを書く方法
これをやってみたい場合は、O’Reilly社から出版されている『Head First Excel』が参考になります。
なぜスプレッドシートでのデータ分析を教えないのでしょうか?
スプレッドシートでのデータ分析には、通常、多大な手作業を要します。パラメータを変更したり、新しいデータセットで分析を実行したりしたい場合、通常はすべてを手作業でやり直さなければなりません。(マクロを作成できることは承知していますが、次の項目をご覧ください。)
また、作業に戻ったり、誰かに分析の詳細を求められたりした際に、スプレッドシートソフトで行った統計解析やグラフ作成の経緯を追跡したり再現したりすることは困難です。
数多くのスプレッドシートソフトが利用可能です。ほとんどの受講者が主要なスプレッドシートソフトとしてExcelを使用しているため、本レッスンではExcelの例を用います。また、無料で利用できるスプレッドシートソフトとしてLibreOfficeもあります。プログラムによってコマンドが多少異なる場合がありますが、基本的な考え方は同じです。
スプレッドシートソフトには、研究者として必要となる機能の多くが盛り込まれています。以下のような用途に活用できます。
- データ入力
- データの整理
- データの抽出と並べ替え
- 統計学
- プロット
スプレッドシートソフトでは、表を用いてデータを表現・表示します。表形式でフォーマットされたデータは本章の主要なテーマでもあり、後続の解析を効率的に行うために、データを標準的な方法で表にまとめる方法について解説します。
課題:隣の人と以下の点について話し合ってください。
- 研究や授業、あるいは自宅でスプレッドシートを使ったことはありますか?
- スプレッドシートではどのような操作を行いますか?
- スプレッドシートはどのような用途に向いていると思いますか?
- スプレッドシートソフトで、うっかり何かをしてイライラしたり、がっかりしたりしたことはありますか?
スプレッドシートの問題点
スプレッドシートはデータ入力に適していますが、実際にはデータ入力以外にも幅広く利用されています。論文用のデータテーブルの作成や、要約統計量の算出、図表の作成などに活用されています。
出版用の表をスプレッドシートで作成するのは最適ではありません。出版用にデータ表をフォーマットする際、主要な要約統計量を本来データとして読み取られることを想定しない形式で提示することが多く、また特殊な書式設定(セルの結合、罫線の設定、見栄えの調整など)を伴うことがよくあります。こうした作業は、文書編集ソフト内で行うことをおすすめします。
後者の2つの用途、すなわち統計量や図表の作成については、注意が必要です。スプレッドシートプログラムはグラフィカルなドラッグ&ドロップ方式であるため、特に統計量や図表の作成により複雑な計算が必要な場合、自身の操作手順を再現することは極めて困難、あるいは不可能に近いと言えます。さらに、スプレッドシートで計算を行う際は、隣接する複数のセルに誤ってわずかに異なる数式を適用してしまいがちです。RやSASのようなコマンドラインベースの統計プログラムを使用する場合、意図的に行わない限り、データセット内のある観測値には計算を適用し、別の観測値には適用しないということは事実上不可能です。
スプレッドシートにおけるデータテーブルの書式設定
質問
- データを効果的に活用するために、スプレッドシートではどのようにデータを整形すればよいでしょうか。
目的
スプレッドシートにおけるデータ入力と書式設定のベストプラクティスについて説明してください。
スプレッドシートで変数と観測値を整理する際は、ベストプラクティスを適用してください。
要点
生データは絶対に変更しないでください。変更を加える前には、必ずコピーを作成してください。
データのクリーニングで行ったすべての手順を、プレーンテキストファイルに記録しておいてください。
整然データの原則に従ってデータを整理してください。
最もよくある間違いは、スプレッドシートソフトを実験ノートのように扱ってしまうことです。つまり、文脈や余白のメモ、データやフィールドの配置に頼って情報を伝えようとすることです。人間であれば(通常は)こうした情報を解釈できますが、コンピュータは同じように情報を捉えません。一つひとつの要素が何を意味するのかをコンピュータに説明しない限り(これは難しい作業です!)、コンピュータはデータ同士がどのように関連し合っているのかを理解することができません。
コンピュータの力を活用すれば、はるかに効率的かつ迅速にデータを管理・分析できます。しかし、その力を活かすには、コンピュータが理解できる形でデータを準備する必要があります(コンピュータは非常に文字通りに解釈するためです)。
だからこそ、最初の予備実験でデータ入力を始める前から、あらかじめ適切にフォーマットされた表を用意しておくことが極めて重要なのです。データの整理は研究プロジェクトの基盤となります。分析を進める過程でデータの扱いが容易になることもあれば、逆に困難になることもあるため、データ入力や実験のセットアップを行う段階で検討しておく価値があります。スプレッドシートではさまざまな方法で設定できますが、その選択によっては、他のプログラムでデータを扱ったり、半年後の自分や共同研究者がデータを扱う際に支障をきたしたりする可能性があります。
注: データ入力やデータ分析に最適なレイアウトやフォーマット(ならびにソフトウェアやインターフェース)は、それぞれ異なる場合があります。この点に留意し、理想的には一方から他方への変換を自動化することが重要です。
解析の記録管理
スプレッドシートを用いてデータの整理や分析を行っていると、元の状態とは大きく異なるスプレッドシートになってしまうことがよくあります。分析結果を再現したり、査読者や指導教員から別の分析を求められた際に何を行ったかを確認したりするためには、
クリーニングまたは解析済みのデータを用いて、新しいファイルを作成してください。元のデータセットは変更しないでください。さもないと、どこから始めたのか分からなくなってしまいます!
クリーニングや解析で行った手順を記録しておいてください。これらの手順は、実験における他の手順と同様に記録してください。データファイルと同じフォルダに保存したプレーンテキストファイルに記録することをお勧めします。
これはスプレッドシートの構成例かもしれません。

これらの原則を、本日の演習で早速実践してみましょう。
本コースではバージョン管理は範囲外となりますが、データのバージョン管理方法については、The Carpentriesの‘Git’に関するレッスンをご参照ください。また、簡単なチュートリアルについてはこちらのブログ記事を、より研究向けの活用事例については Perez-Riverol et al. (2016) をご覧ください。
スプレッドシートでのデータ構造化
データ処理にスプレッドシートソフトを使用する際の基本原則:
- すべての変数を列に配置してください。「weight」や「temperature」のように、測定対象となるものです。
- 各観測値をそれぞれ別の行に記述してください。
- 1つのセルに複数の情報をまとめて入力しないでください。一見すると1つの項目のように見えても、データの利用や並べ替えにおいて、その方法しか選択肢がないのかどうかを考えてみてください。
- 生データはそのままにしておいてください。変更しないでください!
- クリーニング済みのデータをCSV(カンマ区切り値)形式などのテキストベースの形式でエクスポートしてください。これにより、誰でもデータを利用できるようになります。また、ほとんどのデータリポジトリでこの形式が求められています。
例えば、ベルギーのブリュッセルにある複数の病院を受診した患者のデータがあります。そこには、受診日、病院、患者の性別、体重、および血液型が記録されています。
このようにデータを管理したとしたら:

問題は、ABO群とRhesus群が同じBlood型の列にあることです。そのため、A群のすべての観測値を確認したり、ABO群ごとの体重分布を調べたりしたい場合、このデータ構成では困難になります。代わりにABO群とRhesus群を別々の列に配置すれば、はるかに簡単になることがわかります。

データシートを作成する際の重要なルールは、列を変数に、行を観測値に用いることです。
- 列は変数です。
- 行は観測値です。
- セルは個々の値です。
課題:乱雑なデータセットを取り上げ、どのようにクリーニングするかを説明します。
こちらのリンクをクリックして、乱雑なデータセットをダウンロードしてください。
スプレッドシートソフトでデータを開いてください。
2つのタブがあることがわかります。データには、2020年の新型コロナウイルス感染症の第1波および第2波の際に、ブリュッセルのさまざまな病院で記録された多様な臨床変数が含まれています。ご覧の通り、3月の波と11月の波では、データの記録方法が異なります。さて、あなたはこのプロジェクトの担当者となり、データの分析に取りかかりたいと考えています。
隣の人と一緒に、このスプレッドシートの問題点を指摘してください。また、第1波と第2波のタブを整理し、それらを一つのスプレッドシートにまとめるためにどのような手順を踏む必要があるかについても話し合ってください。
重要: 最初のアドバイスを忘れないでください。クリーニング済みデータ用の新しいファイル(またはタブ)を作成し、元の(生の)データは絶対に修正しないでください。
この演習を終えた後、このデータのどこに問題があったのか、そしてどのように修正すべきかについて、グループで話し合います。
課題:データを整理したら、以下の質問に答えてください。
- この研究には男女それぞれ何人が参加しましたか?
- A型、AB型、B型はそれぞれ何例検査されましたか?
- 上記と同様ですが、汚染サンプルを除外した場合はどうなりますか?
- Rhesus +と-はそれぞれ何例検査されましたか?
- これまでに何人の万能ドナー(O-)が検査されましたか?
- AB型男性の平均体重はどれくらいですか?
- 各病院でいくつの検体が検査されましたか?
特にRのスクリプトに関しては、Tidy Dataの論文 Wickham (2014). が優れた参考資料となります。
スプレッドシートでよくある間違い
質問
- スプレッドシートでデータをフォーマットする際によくある課題にはどのようなものがあり、それらをどのように回避できるでしょうか。
目的
- スプレッドシートでよくある書式の問題を認識し、解決します。
要点
- 1つのスプレッドシート内で複数の表を使用することは避けてください。
- データを複数のタブに分散させないようにしてください。
- 0は0として記録してください。
- 欠損データの記録には、適切なヌル値を使用してください。
- 情報を伝えたり、スプレッドシートの見栄えを良くしたりするために、書式設定を使用しないでください。
- コメントは別の列に記入してください。
- 列見出しに記録単位を記載してください。
- セルには情報を1つだけ入力してください。
- 列見出しにスペース、数字、特殊文字は使用しないでください。
- データに特殊文字を含めないようにしてください。
- メタデータは別のプレーンテキストファイルに記録してください。
ご自身のデータだけでなく、共同研究者やインターネットから入手したデータにおいても、注意すべき潜在的なエラーがいくつかあります。こうしたエラーや、その後のデータ分析および結果の解釈に及ぼし得る悪影響を認識していれば、ご自身やプロジェクトメンバーがそれらを回避しようという意欲につながるでしょう。スプレッドシートでのデータ書式設定にわずかな変更を加えるだけでも、データのクリーニングや分析における効率性と信頼性に大きな影響を与える可能性があります。
- 複数の表の使用
- 複数のタブの使用
- ゼロを補完しない
- 問題のある欠損値の使用
- 書式を用いて情報を伝える
- 書式設定でデータシートをきれいに見せる
- セルへのコメントや単位の入力
- セルに複数の情報を入力する
- 問題のあるフィールド名の使用
- データにおける特殊文字の使用
- データテーブルへのメタデータの組み込み
複数の表の使用
よくある手法として、1つのスプレッドシート内に複数のデータテーブルを作成することが挙げられます。これはコンピュータを混乱させるため、絶対に行わないでください!1つのスプレッドシート内に複数のテーブルを作成すると、コンピュータは各行を観測値として扱うため、データ間に誤った関連性を導き出してしまいます。また、同じフィールド名を複数箇所で使用する可能性もあり、そうなるとデータを使いやすい形式に整理することがより困難になります。以下の例がその問題を示しています:

上記の例では、コンピュータは(例えば)4行目を認識し、列A~AFがすべて同一のサンプルを指していると判断します。この行は実際には、4つの異なる採取日(5月29日、6月12日、6月19日、6月26日)それぞれのサンプル1という4つの個別のサンプルと、そのうち2つのサンプルに関する算出済みの要約統計量(平均値(avr)および測定標準誤差(SEM))を表しています。他の行も同様に問題があります。
複数のタブの使用
では、ワークブックのタブはどうでしょうか。データを整理する簡単な方法に思えますよね。ですが、一概にそうとは言えません。タブを追加してしまうと、データに存在する関係性をコンピュータに認識させることができなくなります(この関係性を確保するには、スプレッドシートアプリケーション固有の関数やスクリプトを導入する必要があります)。例えば、測定を行う日ごとに個別のタブを作成したとしましょう。
これは2つの理由から、好ましくない手法です。
測定を行うたびに新しいタブでデータの記録を始めると、誤ってデータに不整合が生じやすくなります。また、
たとえあらゆる不整合の混入を防げたとしても、これらのデータを一つのデータテーブルに統合する必要があるため、分析前に余計な手間がかかることになります。タブをどのように統合するかをコンピュータに明示的に指示しなければならず、もしタブの書式が統一されていない場合は、手作業で行わなければならないことさえあります。
次にデータを入力する際、別のタブやテーブルを作成しようとしたら、元のスプレッドシートに列を一つ追加するだけで、このタブを省けないか考えてみてください。乱雑なデータファイルの例では複数のタブを使用しましたが、これでデータを再構成し、タブを統合する方法がお分かりいただけたかと思います。
実験が進むにつれて、データシートが非常に長くなることがあります。スプレッドシートの上部に見出しが表示されないと、データの入力が困難になります。ただし、見出し行を繰り返さないでください。見出し行がデータに紛れ込みやすく、後々問題の原因となることがあります。代わりに、列見出しを固定すれば、行数の多いスプレッドシートでも見出しを表示したままにできます。
ゼロを補完しない
何かを測定する場合、通常はゼロになることがあります。例えば、調査でウサギが観察された回数などです。ほとんどがゼロなのに、わざわざその列に「0」と書き込む必要があるでしょうか。
ただし、スプレッドシートにおける「0」と「空白セル」は異なります。コンピュータにとって、0は実際にはデータです。測定またはカウントされたものです。一方、空白セルは測定されなかったことを意味し、コンピュータはそれを未知の値(ヌル値または欠損値とも呼ばれます)として解釈します。
スプレッドシートや統計ソフトでは、ゼロとして意図した空白セルが誤って解釈される可能性があります。観測値を入力しないことで、そのデータを「不明」または「欠損(null)」として扱うようコンピュータに指示することになります。これにより、後続の計算や分析で問題が生じることがあります。例えば、1つの欠損値を含む数値群の平均値は、常に欠損値となります(コンピュータが欠損した観測値の値を推測できないためです)。このため、ゼロはゼロとして、真に欠損しているデータは欠損値として記録することが極めて重要です。
問題のある欠損値の使用
例: -999やその他の数値(または0)を用いて欠損データを表す方法。
解答:
データセット内で欠損値がそれぞれ異なる形式で表現されるのには、いくつかの理由があります。測定機器から紛らわしい欠損値が自動的に記録される場合があります。その場合はどうすることもできませんが、分析前にOpenRefineのようなツールを用いてデータクリーニングを行うことで対処可能です。また、データが存在しない理由に応じて、異なる欠損値が使い分けられることもあります。これは把握しておくべき重要な情報ですが、実質的に1つの列で2つの情報を表現していることになります。書式を用いて情報を伝える場合と同様に、ここでも「data_missing」のような新しい列を作成し、その列で異なる理由を記録するのが良いでしょう。
理由はどうあれ、未知のデータや欠損データが-999、999、または0として記録されているのは問題です。
多くの統計ソフトでは、これらが欠損値(null値)を表す意図であるとは認識されません。これらの値がどのように解釈されるかは、データの解析に使用するソフトウェアによって異なります。明確に定義され、一貫したnull指標を用いることが不可欠です。
空白(ほとんどのアプリケーション)やNA(Rの場合)が適しています。 White et al. (2013) は記事の中で、さまざまなソフトウェアアプリケーションにおけるヌル値の表記として適切な選択肢について解説しています。

書式を用いて情報を伝える
例: 解析から除外すべきセル、行、または列をハイライトしたり、データの区切りを示すために空行を設けたりすることです。

解決策: 除外すべきデータをエンコードするための新しいフィールドを作成します。

書式設定を用いてデータシートの見栄えを整える
例: セルの結合。
解決策: 注意しないと、ワークシートの見栄えを良くするために書式設定を行うと、コンピュータがデータ内の関連性を認識できなくなる恐れがあります。セルを結合すると、統計ソフトでデータを読み取れなくなります。データを整理する際にセルを結合する必要がないよう、データの再構成を検討してください。
セルへのコメントや単位の入力
ほとんどの解析ソフトウェアはExcelやLibreOfficeのコメントを認識できず、データセル内に記述されたコメントを誤認してしまいます。書式設定について前述した通り、セルにメモを追加する必要がある場合は、別のフィールドを作成してください。同様に、セル内に単位を含めないでください。理想的には、1つの列に配置するすべての測定値は同じ単位であるべきですが、何らかの理由で異なる場合は、別のフィールドを作成し、そのセルの単位を指定してください。
セルに複数の情報を入力する
例: ABO型とRh型のグループを、A+、B+、A-…のように1つのセルに記録します。
解決策: 1つのセルに複数の情報を入れないでください。1つのセルに複数の情報を入れると、データの分析方法が制限されてしまいます。これら両方の測定値が必要な場合は、それらの情報を含めるようにデータシートを設計してください。例えば、ABO群用とRh群用にそれぞれ1列ずつ設けてください。
問題のあるフィールド名の使用
説明的なフィールド名を選択してください。ただし、スペースや数字、いかなる種類の特殊文字も含まないようご注意ください。スペースは、空白を区切り文字として使用するパーサーによって誤認される可能性があります。また、一部のプログラムでは、数字で始まるテキスト文字列のフィールド名を処理できない場合があります。
アンダースコア(_)はスペースに代わる有効な選択肢です。可読性を高めるため、名称はキャメルケース(例:ExampleFileName)で記述することを検討してください。現時点では意味が通じる略語であっても、6か月後には必ずしもそうとは限らない点に留意してください。ただし、名称を過度に長くしすぎないようにしましょう。フィールド名に単位を含めることで混乱を防ぎ、他者がフィールドの内容を容易に解釈できるようになります。
例
| 適切な名称 | 優れた代替案 | 避ける |
|---|---|---|
| Max_temp_C | MaxTemp | Maximum Temp (°C) |
| Precipitation_mm | Precipitation | precmm |
| Mean_year_growth | MeanYearGrowth | Mean growth/year |
| sex | sex | M/F |
| weight | weight | w. |
| cell_type | CellType | Cell Type |
| Observation_01 | first_observation | 1st Obs |
データにおける特殊文字の使用
例: メモを作成する際、スプレッドシートソフトをワードプロセッサのように扱い、例えばWordや他のアプリケーションからデータを直接コピーすることがあります。
解決策: これはよくある手法です。例えば、セルに長い文章を記述する際、スプレッドシートに改行やエムダッシュなどを含めることがよくあります。また、Wordなどのアプリケーションからデータをコピーすると、書式や特殊な非標準文字(左揃え・右揃えの引用符など)が含まれてしまいます。こうしたデータをコーディング・統計環境やリレーショナルデータベースにエクスポートすると、行が途中で切れたり、エンコーディングエラーが発生したりといった問題が生じる可能性があります。
一般的なベストプラクティスとしては、改行やタブ、垂直タブなどの文字を追加しないことが推奨されます。言い換えれば、テキストセルはテキストとスペースしか入力できないシンプルなウェブフォームであるかのように扱ってください。
データテーブルへのメタデータの組み込み
例: データテーブルの上部または下部に、列の意味や単位、例外事項などを説明する凡例を追加します。
解決策: データに関する情報(「メタデータ」)を記録しておくことは不可欠です。データの収集や分析を行っている最中はデータセットに精通しているかもしれませんが、数ヶ月、1年、あるいはそれ以上の時間が経過した後に、例えば変数「sglmemgp」がグループの単一メンバーを意味することや、変数の変換や派生変数の作成に使用した正確なアルゴリズムを、依然として覚えていられる可能性は低いでしょう。
また、他者があなたのデータを精査したり利用したりしたいと考える理由は数多くあります。得られた知見を理解するため、その妥当性を検証するため、投稿した論文を査読するため、結果を再現するため、同様の研究を設計するため、あるいは他者がアクセスし再利用できるようデータをアーカイブするためなどです。デジタルデータは定義上、機械可読ですが、その意味を理解するのは人間の役割です。研究の収集・分析段階においてデータを文書化することの重要性は、とりわけその研究が学術的記録の一部となる場合には、いくら強調してもしすぎることはありません。
ただし、メタデータはデータファイル自体に含めるべきではありません。論文の表や補足ファイルとは異なり、メタデータ(凡例形式)はデータではないため、データファイルに含めるべきではありません。また、これを含めると、コンピュータプログラムによるデータファイルの解釈に支障をきたす可能性があります。むしろ、メタデータはデータファイルと同じディレクトリに別ファイルとして保存し、できればデータファイルとの関連性が明確にわかる名称を付けたプレーンテキスト形式で保存してください。メタデータファイルは自由形式であるため、文書化は重要でありながらデータファイルのフォーマットを乱す可能性のあるコメントや単位、欠損値の符号化方法などの情報を記述することも可能です。
さらに、ファイルまたはデータベースレベルのメタデータでは、データセットを構成するファイル同士の関係性やフォーマット、およびそれらが以前のファイルに優先するか、あるいはそれらに取って代わられるかが記述されます。フォルダレベルの readme.txt ファイルは、プロジェクト内のすべてのファイルとフォルダを記録する定番の手法です。
(メタデータのテキストは、エディンバラ大学のEDINAおよびData Libraryによるオンラインコース「Research Data MANTRA」を基に作成されています。MANTRAはクリエイティブ・コモンズ 表示 4.0 国際ライセンスの下でライセンスされています。)
データのエクスポート
質問
- スプレッドシートから後続のアプリケーションで活用しやすい形でデータをエクスポートするにはどうすればよいでしょうか。
目的
- スプレッドシートのデータを汎用的なファイル形式で保存します。
- スプレッドシートからCSVファイルにデータをエクスポートします。
要点
一般的なスプレッドシート形式で保存されたデータは、データ分析ソフトウェアで正しく読み込まれないことが多く、その結果、データにエラーが生じることがあります。
スプレッドシートからCSVやTSVなどの形式にデータをエクスポートすることで、ほとんどのプログラムで一貫して利用できる形式になります。
分析に使用するデータを、Excelのデフォルトファイル形式(Excelのバージョンに応じて*.xlsまたは*.xlsx)で保存するのは得策ではありません。なぜでしょうか?
これはプロプライエタリな形式であり、将来的には、そのファイルを開くための技術が存在しなくなる(または非常にまれになる)ことで、ファイルを開くことが不可能になったり、少なくとも不便になったりする可能性があるためです。
他のスプレッドシートソフトでは、独自のExcel形式で保存されたファイルを開けない場合があります。
Excelのバージョンによってデータの扱いが異なる場合があり、それが不整合を招くことがあります。日付は、データ保存における不整合のよく知られた例です。
最後に、データリポジトリへのデータ寄託を求める学術誌や助成機関が増えており、その多くはExcel形式を受け付けていません。以下で説明する形式のいずれかである必要があります。
上記の点は、LibreOffice / Open Officeで使用されるオープンデータ形式などの他の形式にも当てはまります。これらの形式は静的なものではなく、ソフトウェアパッケージによって解析方法が異なります。
汎用的でオープン、かつ静的な形式でデータを保存することは、この問題の解決に役立ちます。タブ区切り(タブ区切り値またはTSV)やカンマ区切り(カンマ区切り値またはCSV)を試してみてください。CSVファイルは、列がカンマで区切られたプレーンテキストファイルであり、そのため「カンマ区切り値」またはCSVと呼ばれます。ExcelやSPSSなどのファイルに対するCSVファイルの利点は、TextEditやメモ帳のようなプレーンテキストエディタを含め、ほぼあらゆるソフトウェアでCSVファイルを開いて読み取れる点にあります。また、CSVファイル内のデータは、SQLiteやRなどの他の形式や環境へも簡単にインポートできます。CSVファイルを扱う際は、特定の高価なプログラムの特定のバージョンに縛られることがないため、最大限の移植性と永続性を確保する上で優れた形式と言えます。ほとんどの表計算ソフトでは、ファイルエクスポート時に警告が表示される場合がありますが、CSVのような区切り文字付きテキスト形式へ簡単に保存できます。
Excelで開いたファイルをCSV形式で保存するには:
- 上部メニューから「ファイル」→「名前を付けて保存」を選択してください。
- 「書式」フィールドで、一覧から「カンマ区切り値」(
*.csv)を選択してください。 - ファイル名と保存先を再確認し、「保存」をクリックしてください。
後方互換性に関する重要な注意点:CSVファイルはExcelで開くことができます!

Rとxlsに関する注記:
xlsファイル(およびGoogleスプレッドシート)を読み込めるRパッケージがあります。xlsドキュメント内の異なるワークシートにアクセスすることさえ可能です。
しかし
- これらの中には、Windowsでしか動作しないものもあります。
- これは、(単純ながらも手作業による)
csvへのエクスポートを、データ分析用のRコードにおけるさらなる複雑さや依存関係に置き換えることに相当します。 - データフォーマットのベストプラクティスは引き続き適用されます。
-
csv(あるいはそれに類するもの)では不十分な、本当に納得のいく理由があるのでしょうか?
カンマに関する注意点
一部のデータセットでは、データ値自体にカンマ(,)が含まれている場合があります。その場合、使用するソフトウェア(Excelを含む)では、データが列単位で誤って表示される可能性が高くなります。これは、データ値の一部であるカンマが区切り文字として解釈されるためです。
例えば、データは次のようになっているかもしれません。
species_id,genus,species,taxa
AB,Amphispiza,bilineata,Bird
AH,Ammospermophilus,harrisi,Rodent, not censused
AS,Ammodramus,savannarum,Bird
BA,Baiomys,taylori,Rodent
レコード
AH,Ammospermophilus,harrisi,Rodent, not censused の
taxa の値にはカンマ (Rodent, not censused)
が含まれています。上記の内容を
Excel(または他のスプレッドシートソフト)で読み込もうとすると、次のような結果になります。

taxa の値は(1つの列 D
にまとめられるのではなく)2つの列に分割されました。これにより、さらに多くのエラーが連鎖的に発生する可能性があります。例えば、追加された列は多くの欠損値を含む列(かつ適切なヘッダーのない列)として解釈されます。加えて、3行目のレコード(つまり
‘taxa’ の値にカンマが含まれていた行)における列 D
の値が、現在は正しくありません。
データをcsv形式で保存したい場合で、データ値にカンマが含まれる可能性があるときは、値を引用符(““)で囲むことで、前述の問題を回避できます。このルールを適用すると、データは次のようになります。
species_id,genus,species,taxa
"AB","Amphispiza","bilineata","Bird"
"AH","Ammospermophilus","harrisi","Rodent, not censused"
"AS","Ammodramus","savannarum","Bird"
"BA","Baiomys","taylori","Rodent"
このファイルをExcelでcsvとして開いても、余分な列が追加されることはありません。Excelでは、引用符の外側にあるカンマのみを区切り文字として使用するためです。
あるいは、カンマを含むデータを扱う場合は、スプレッドシートで作業する際に別の区切り文字を使用する必要があるでしょう1。この場合は、区切り文字としてタブを使用し、TSVファイルを扱うことを検討してください。TSVファイルは、CSVファイルと同様の方法でスプレッドシートソフトからエクスポートできます。
データ値が「““」で囲まれておらず、区切り文字とデータ値の一部の両方にカンマが使用されている既存のデータセットを扱っている場合、データクリーニングにおいて重大な問題に直面している可能性があります。扱うデータセットに数百から数千件のレコードが含まれている場合、手動でクリーニングを行う(データ値からカンマを削除するか、値を引用符”” で囲む)作業は、膨大な時間を要するだけでなく、誤って多くのエラーを混入させてしまう恐れもあります。
データセットのクリーニングは、多くの科学分野における主要な課題の一つです。その手法は、ほぼ常に具体的な状況によって異なります。しかし、例えばスクリプトを記述して実行するなど、自動化された方法でデータをクリーニングすることは良い習慣です。PythonとRのレッスンでは、適切なスクリプトを作成するためのスキルを身につける基礎を学ぶことができます。
概要

典型的なデータ分析のワークフローを上の図に示しています。ここでは、データの変換、可視化、モデリングが繰り返し行われます。この反復作業は、データが理解されるまで何度も繰り返されます。しかし実際の現場では、実際に分析して理解することよりも、データのクリーニングや準備に大半の時間が費やされることが多いのが現状です。
変換・可視化・モデリングのサイクルを迅速に数回繰り返すアジャイルなデータ分析ワークフローが実現可能なのは、データが予測可能な形式で整備されており、データを直接確認したり修正したりすることなく内容を推論できる場合に限られます。
- 適切なデータ整理は、あらゆる研究プロジェクトの基盤となります。
これは、カンマを小数点区切り文字として使用する欧州諸国において特に当てはまります。このような場合、CSVファイルのデフォルトの値区切り文字はセミコロン(;)になるか、あるいは値が一律に引用符で囲まれます。↩︎
Content from RとRStudio
最終更新日:2026-11-07 | ページの編集
所要時間: 30分
概要
質問
- RとRStudioとは何ですか?
目的
- RStudioのスクリプト、コンソール、環境、プロットの各ペインの目的について説明してください。
- 一連の解析のためのファイルとディレクトリをRプロジェクトとして整理し、作業ディレクトリの目的を理解します。
- Rの関数に関する詳細情報を調べるには、RStudio内蔵のヘルプインターフェースをご利用ください。
- Rユーザーコミュニティでトラブルシューティングを行う際に、十分な情報を提供する方法を説明します。
このエピソードは、Data Carpentriesの「生態学者のためのRによるデータ分析と可視化」レッスンに基づいています。
Rとは何ですか?RStudioとは何ですか?
Rという用語は、プログラミング言語、統計計算環境、およびそれを用いて記述されたスクリプトを解釈するソフトウェアを指して用いられます。
RStudioは現在、Rスクリプトを記述するだけでなく、Rソフトウェアを操作するためにも非常に広く利用されているツールです1。正常に動作させるにはRが必要であるため、お使いのコンピュータにRとRStudioの両方がインストールされている必要があります。
RStudio IDEチートシートには、ここで扱う内容よりもはるかに多くの情報が記載されていますが、キーボードショートカットを覚えたり、新機能を発見したりするのに役立ちます。
なぜRを学ぶのでしょうか?
Rはあちこちを指し示したりクリックしたりすることがあまりありません。それは良いことです。
他のソフトウェアに比べて習得までのハードルは高いかもしれませんが、Rでは分析結果が一連のポインティング&クリック操作の記憶に依存せず、一連の記述式コマンドに基づいて得られるという点は大きなメリットです。そのため、さらにデータを収集して分析をやり直したい場合でも、結果を得るためにどのボタンをどの順番でクリックしたかを思い出す必要はありません。単にスクリプトを再実行すればよいのです。
スクリプトを活用することで、分析で用いた手順が明確になります。また、記述したコードを他の人が確認し、フィードバックをくれたりミスを見つけたりすることも可能です。
スクリプトを活用することで、自身の作業に対する理解を深めることができ、使用する手法の習得や理解も容易になります。
Rコードは再現性に優れています。
再現性とは、同じ解析コードを用いれば、他の誰か(将来の自分を含む)が同じデータセットから同じ結果を得られることを意味します。
Rは他のツールと連携し、コードから原稿やレポートを作成できます。さらにデータを収集したり、データセットの誤りを修正したりすると、原稿やレポート内の図や統計検定が自動的に更新されます。
再現性のある解析を求める学術誌や助成機関が増えており、Rを習得していれば、こうした要件への対応において優位に立つことができます。
Rは学際的で拡張性に優れています。
機能を拡張するためにインストール可能な10,000以上のパッケージ2を備えたRは、多くの科学分野の統計的手法を組み合わせ、データ解析に必要な分析フレームワークに最適に適合させることができる枠組みを提供しています。例えば、Rには画像解析、GIS、時系列解析、集団遺伝学など、多岐にわたるパッケージが用意されています。

Rはあらゆる形式や規模のデータを扱えます。
Rで習得するスキルは、データセットの規模に応じて容易に応用できます。データセットが数百行であっても数百万行であっても、それほど大きな違いはありません。
Rはデータ分析向けに設計されています。欠損データや統計的因子の取り扱いを容易にする、専用のデータ構造やデータ型を備えています。
Rは、お使いのコンピュータ上やウェブ上のスプレッドシート、データベース、その他多くのデータ形式に接続できます。
Rには大規模で親しみやすいコミュニティがあります。
毎日、数千人もの人々がRを利用しています。その多くは、メーリングリストやStack Overflowなどのウェブサイト、あるいはRStudioコミュニティを通じて、快くサポートしてくれます。こうした広範なユーザーコミュニティは、バイオインフォマティクスなどの専門分野にも及んでいます。Rコミュニティのこうした一部を構成するのがBioconductorであり、これは「現在および新興の生物学的アッセイによるデータ」の解析と理解を目的とした科学プロジェクトです。本ワークショップはBioconductorコミュニティのメンバーによって開発されました。Bioconductorの詳細については、関連ワークショップである“The Bioconductor Project”をご覧ください。
RStudioの基本操作
まずは、Rを扱うための統合開発環境(IDE)であるRStudioについて学んでいきましょう。
オープンソース製品であるRStudio IDEは、Affero General Public License (AGPL) v3に基づき無料で提供されています。また、RStudio IDEはPosit, Inc.から商用ライセンスおよび優先メールサポート付きで提供されています。
コードの記述、コンピュータ上のファイルの操作、作成予定の変数の確認、および生成するプロットの可視化には、RStudio IDEを使用します。RStudioは、本ワークショップでは扱わない他の用途(例:バージョン管理、パッケージ開発、Shinyアプリの作成)にも利用できます。

RStudioのウィンドウは4つの「ペイン」に分かれています。
- スクリプトやドキュメントのソース(デフォルトレイアウトでは左上)
- 環境/履歴(右上)、
- Files/Plots/Packages/Help/Viewer(右下)、および
- Rのコンソール(左下)。
これらのペインの配置や内容はカスタマイズ可能です(メニューのTools -> Global Options -> Pane Layoutを参照してください)。
RStudioを利用するメリットの一つは、コードを記述するために必要なすべての情報が一つのウィンドウで確認できる点です。さらに、Rでの開発で使用する主要なファイル形式に対する多くのショートカット、自動補完、ハイライト機能により、RStudioを使えばより簡単に、かつミスの少ない入力が可能になります。
セットアップ
関連するデータ、解析、テキストを、作業ディレクトリと呼ばれる単一のフォルダにまとめて管理するのが一般的です。このフォルダ内のすべてのスクリプトでは、プロジェクト内のどこにファイルがあるかを示す相対パスを使用できます(特定のコンピュータ上の場所を指す絶対パスとは異なります)。このように管理することで、基盤となるスクリプトが引き続き動作するかどうかを心配することなく、コンピュータ上でプロジェクトを移動させたり、他者と共有したりすることが格段に容易になります。
RStudioでは、「プロジェクト」インターフェースを通じてこれを実現するための便利なツール群が提供されています。この機能は、作業ディレクトリを作成するだけでなく、その場所を記憶し(これにより素早く移動できます)、オプションでカスタム設定や開いているファイルを保持することで、中断後の作業再開を容易にします。以下の手順に従って、このチュートリアル用の「Rプロジェクト」を作成してください。
- RStudioを起動します。
-
Fileメニューから、New projectをクリックしてください。New directoryを選択し、次にNew projectを選択してください。 - この新しいフォルダ(または「ディレクトリ」)に名前を入力し、適切な場所を指定してください。これが今回のセッション(またはコース全体)の作業ディレクトリとなります(例:
bioc-intro)。 -
Create projectをクリックしてください。 - (任意)RStudioの環境設定で、ワークスペースを保存しないよう「Never」に設定してください。
RStudioのデフォルト設定は概ね問題なく動作しますが、特に大規模なデータセットを扱う場合、ワークスペースを.RDataに保存するのは手間がかかることがあります。この機能を無効にするには、「ツール」→「グローバルオプション」に進み、「終了時にワークスペースを.RDataに保存」のオプションで「なし」を選択してください。

Windowsと他のOS間の文字エンコーディングの問題を避けるため、デフォルトでUTF-8を設定します。

作業ディレクトリの整理
プロジェクト間で一貫したフォルダ構成を採用することで、整理された状態を維持でき、将来的にファイルの検索や整理も容易になります。これは、複数のプロジェクトを抱えている場合に特に有用です。一般的には、スクリプト、データ、ドキュメント用のディレクトリ(フォルダ)を作成するとよいでしょう。
-
data/このフォルダを、生データや、特定の解析のために作成する中間データセットの保存先に使用してください。透明性と来歴を確保するため、生データのコピーは常に参照可能な状態で保持し、データのクリーニングや前処理は可能な限りプログラム的に(つまり、手動ではなくスクリプトを用いて)行うようにしてください。また、生データと処理済みデータを分離しておくことも推奨されます。例えば、scripts/01.preprocess.tree_survey.Rスクリプトによって生成されたdata/processed/tree.survey.csvファイルとは別に、data/raw/tree_survey.plot1.txtや...plot2.txtといったファイルを保持するといった方法が考えられます。 -
documents/ここには、アウトラインや下書き、その他のテキストを保存します。 -
scripts/(またはsrc) ここには、各種解析やプロット用のRスクリプトを保存し、場合によっては関数用の別フォルダを設けることになります(これについては後述します)。
プロジェクトのニーズに応じて、追加のディレクトリやサブディレクトリが必要になる場合もありますが、これらが作業ディレクトリの骨格となるべきです。

本コースでは、生データを保存するためのdata/フォルダが必要になります。また、データをCSVファイルとしてエクスポートする方法を学ぶ際にはdata_output/を、保存する図にはfig_output/フォルダを使用します。
課題:プロジェクトのディレクトリ構成を作成してください。
画面右側のFilesタブでNew Folderをクリックし、新しく作成した作業ディレクトリ(例:~/bioc-intro/data)内にdataという名前のフォルダを作成してください。(または、Rコンソールでdir.create("data")と入力してください。)これらの操作を繰り返し、data_output/とfig_outputのフォルダを作成してください。
使用するファイルは1つだけですし、その方が手間がかからないため、スクリプトは作業ディレクトリのルートに配置しておきます。
作業ディレクトリは、現在以下のようになっているはずです。

プロジェクト管理は、もちろんバイオインフォマティクスのプロジェクトにも適用可能です3。William Noble ( Noble (2009) ) は、次のようなディレクトリ構成を提案しています。
ディレクトリ名は大きめの書体で、ファイル名は小さめの書体で表記されています。ここではファイルの一部のみが表示されています。なお、日付は時系列順に並べ替えられるよう
<year>-<month>-<day>形式で表記されています。ソースコードsrc/ms-analysis.cはコンパイルされてbin/ms-analysisが生成され、doc/ms-analysis.htmlにドキュメントが記載されています。データディレクトリ内のREADMEファイルには、誰が、どのURLから、どの日付にデータファイルをダウンロードしたかが指定されています。ドライバスクリプトresults/2009-01-15/runallは、3つの交差検証分割に対応する3つのサブディレクトリ(split1、split2、split3)を自動的に生成します。bin/parse-sqt.pyスクリプトは、両方のrunallドライバスクリプトから呼び出されます。

適切に定義され、十分に文書化されたプロジェクトディレクトリにおいて最も重要な点は、プロジェクトに不慣れな人4が
プロジェクトの内容、利用可能なデータ、実施された解析、および得られた結果を理解し、そして何よりも重要なのが
分析を再度実施します ― 新しいデータを用いるか、あるいは分析パラメータの一部を変更して行います。
作業ディレクトリ
作業ディレクトリは、理解しておくべき重要な概念です。これは、Rがファイルを探し、保存する場所のことです。プロジェクトのコードを記述する際は、作業ディレクトリのルートを基準としてファイルを参照し、この構造内にあるファイルのみを使用するようにしてください。
RStudioのプロジェクトを利用すれば、この操作は簡単になり、作業ディレクトリも適切に設定されます。確認が必要な場合は、getwd()を使用できます。何らかの理由で作業ディレクトリが意図した状態になっていない場合は、RStudioのインターフェースで、作業ディレクトリがあるべき場所のファイルブラウザに移動し、青い歯車アイコンMoreをクリックしてSet As Working Directoryを選択することで、変更できます。あるいは、setwd("/path/to/working/directory")を使用して作業ディレクトリをリセットすることも可能です。ただし、この行をスクリプトに含めないでください。他の人のコンピュータでは動作しなくなるためです。
例
以下のスキーマは、dataおよびfig_outputというサブディレクトリを持つ作業ディレクトリbioc-introと、後者にある2つのファイルを表しています。
bioc-intro/data/
/fig_output/fig1.pdf
/fig_output/fig2.png
作業ディレクトリにいる場合は、相対パス
bioc-intro/fig_output/fig1.pdf または絶対パス
/home/user/bioc-intro/fig_output/fig1.pdf を用いて
fig1.pdf ファイルを参照できます。
dataディレクトリにいる場合は、相対パス../fig_output/fig1.pdfまたは同じ絶対パス/home/user/bioc-intro/fig_output/fig1.pdfを使用します。
Rの操作
プログラミングの基本は、コンピュータが従うべき指示を書き出し、その指示に従うようコンピュータに指示することです。Rで指示を書き出す、つまりコードを書くのは、Rがコンピュータと人間の双方が理解できる共通言語だからです。これらの指示をコマンドと呼び、それらのコマンドを実行(ランとも呼びます)することで、コンピュータにその指示に従わせます。
Rを操作する方法には、主にコンソールを使う方法と、スクリプト(コードを記述したプレーンテキストファイル)を使う方法の2通りがあります。コンソールペイン(RStudioでは左下のパネル)は、R言語で記述したコマンドを入力し、コンピュータで即座に実行できる場所です。また、実行したコマンドの結果が表示される場所でもあります。コンソールに直接コマンドを入力してEnterを押すことで実行できますが、セッションを閉じるとそれらの内容は消えてしまいます。
コードとワークフローの再現性を確保するためには、実行したいコマンドをスクリプトエディタに入力し、スクリプトを保存しておくのが得策です。そうすることで、行った操作の完全な記録が残り、誰でも(将来の自分も含めて!)自分のコンピュータで簡単に結果を再現できるようになります。ただし、単にスクリプトにコマンドを入力しただけでは、自動的に実行されるわけではない点に注意してください。実行するには、依然としてコンソールに送信する必要があります。
RStudioでは、Ctrl + Enter
のショートカットを使用して、スクリプトエディタから直接コマンドを実行できます(Macでは
Cmd + Return
でも動作します)。Ctrl + Enter
を押すと、スクリプトの現在の行(カーソルで示されている行)にあるコマンド、または現在選択されているテキスト内のすべてのコマンドがコンソールに送信され、実行されます。その他のキーボードショートカットについては、こちらのRStudio
IDEに関するRStudioチートシートをご覧ください。
解析の途中で、スクリプトに必ずしも記録を残すことなく、変数の中身やオブジェクトの構造を確認したい場合があります。これらのコマンドは、コンソールで直接入力して実行できます。RStudioでは、Ctrl
+ 1 および Ctrl + 2
というショートカットを使用して、スクリプトとコンソールペインの間を切り替えることができます。
Rがコマンドを受け付けられる状態になると、Rコンソールに>プロンプトが表示されます。コマンド(入力、コピー&ペースト、またはスクリプトエディタからCtrl
+
Enterを使用して送信)を受け取ると、Rはその実行を試み、準備が整い次第、結果を表示して新しい>プロンプトに戻り、次のコマンドを待機します。
Rがまだ完了しておらず、さらなるデータの入力を待っている状態の場合、コンソールに+というプロンプトが表示されます。これは、完全なコマンドの入力が完了していないことを意味します。これは、括弧や引用符を「閉じて」いない、つまり左括弧の数と右括弧の数、あるいは引用符の開きと閉じの数が一致していないためです。このような状態で、コマンドの入力が完了したと思った場合は、コンソールウィンドウ内をクリックしてEscを押してください。これにより、未完了のコマンドがキャンセルされ、>プロンプトに戻ります。
コース期間中および終了後に、さらに詳しく学ぶにはどうすればよいでしょうか?
本コースで扱う内容を通じて、ご自身の研究のためにRを用いてデータを分析する方法の基礎に触れていただけます。ただし、データセットのクリーニングや統計手法の活用、あるいは見栄えの良いグラフの作成といった高度な操作を行うには、さらなる学習が必要です5。他のツールと同様に、Rに習熟し効率的に使いこなせるようになる最善の方法は、実際の研究課題に取り組むために活用することです。初心者の方にとって、スクリプトを一から書くのは気が遠くなるように感じられるかもしれません。また、多くの人がコードをオンラインで公開しているため、既存のコードを目的に合わせて修正して利用する方が、よりスムーズに学習を始められるでしょう。

ヘルプの参照
RStudio内蔵のヘルプインターフェースを使用して、R関数に関する詳細情報を検索してください。

ヘルプを得る最も迅速な方法の一つは、RStudioのヘルプインターフェースを利用することです。このパネルはデフォルトでRStudioの右下に配置されています。スクリーンショットにあるように、「Mean」という単語を入力すると、RStudioは関心を持ちそうな候補をいくつか提示しようとします。その後、その説明が表示ウィンドウに表示されます。
使いたい関数の名前は分かっているのですが、使い方が分かりません。
特定の関数(例えば
barplot())についてヘルプが必要な場合は、次のように入力してください。
R
?barplot
引数の名前を再確認したいだけであれば、以下をご利用いただけます。
R
args(lm)
Xという処理を行う関数を使いたいのですが、きっとそういう関数があるはずなのに、どれなのか分かりません……。
特定の処理を行う関数をお探しの場合は、ダブルクエスチョンマーク
?? で呼び出す help.search()
関数をご利用いただけます。ただし、この関数はインストール済みパッケージの中から、検索内容に一致するヘルプページを検索するだけです。
R
??kruskal
お探しの情報が見つからない場合は、利用可能なすべてのパッケージのヘルプファイルを検索できる rdocumentation.org のウェブサイトをご利用ください。
最後に、Googleやインターネットで「R <task>」と汎用的に検索すれば、多くの場合、適切なパッケージのドキュメントや、すでに誰かが同じ質問をしている役立つフォーラムに辿り着くでしょう。
行き詰まってしまいました… 意味のわからないエラーメッセージが表示されます。
まずはエラーメッセージをGoogleで検索してみましょう。ただし、パッケージ開発者はRが提供するエラー処理機能を利用していることが多いため、必ずしもうまく機能するとは限りません。結果として、問題の診断にあまり役立たない一般的なエラーメッセージ(例:「subscript out of bounds」)が表示されることがあります。メッセージが非常に汎用的な場合は、使用している関数やパッケージの名前も検索クエリに含めてみてください。
ただし、Stack
Overflowも確認してみてください。[r]タグを使って検索してください。ほとんどの質問にはすでに回答が付いていますが、検索時に適切な言葉を選んで答えを見つけ出すのがポイントです:
https://stackoverflow.com/questions/tagged/r
R入門は、プログラミング経験が少ない方には内容が難解に感じられることもありますが、R言語の基礎を理解するには最適な教材です。
R FAQは内容が濃密で専門的ですが、有用な情報が満載です。
助けを求める
誰かに助けてもらうための鍵は、相手に問題を迅速に把握してもらうことです。問題がどこにあるのかをできるだけ簡単に特定できるようにしてください。
問題を説明する際は、適切な言葉を使うようにしてください。例えば、「パッケージ」と「ライブラリ」は別物です。ほとんどの人は意図を理解してくれるでしょうが、この意味の違いに強いこだわりを持つ人もいます。重要なのは、助けようとしている人を混乱させてしまう可能性があるということです。問題を説明する際は、可能な限り正確に記述してください。
可能であれば、うまくいかない箇所をシンプルな再現可能な例に絞り込んでください。50,000行10,000列のデータフレームではなく、ごく小規模なデータフレームで問題を再現できる場合は、その小規模なデータと問題の説明を添えてください。適宜、ご自身の分野外の人でも質問内容を理解できるよう、記述を一般化するよう心がけてください。例えば、実際のデータセットのサブセットを使うのではなく、小規模な(3列5行の)汎用的なデータセットを作成してください。再現可能な例の書き方の詳細については、Hadley Wickham氏によるこちらの記事を参照してください。
オブジェクトを他の人と共有する場合、サイズが比較的小さければ、dput()関数を使用できます。これにより、メモリ上のものと全く同じオブジェクトを再現できるRコードが出力されます。
R
## iris is an example data frame that comes with R and head() is a
## function that returns the first part of the data frame
dput(head(iris))
出力
structure(list(Sepal.Length = c(5.1, 4.9, 4.7, 4.6, 5, 5.4),
Sepal.Width = c(3.5, 3, 3.2, 3.1, 3.6, 3.9), Petal.Length = c(1.4,
1.4, 1.3, 1.5, 1.4, 1.7), Petal.Width = c(0.2, 0.2, 0.2,
0.2, 0.2, 0.4), Species = structure(c(1L, 1L, 1L, 1L, 1L,
1L), levels = c("setosa", "versicolor", "virginica"), class = "factor")), row.names = c(NA,
6L), class = "data.frame")
オブジェクトのサイズが大きい場合は、生ファイル(つまりCSVファイル)と、エラーが発生する箇所までのスクリプト(問題に関係のない部分はすべて削除したもの)を併せて提供してください。あるいは、特に質問内容がデータフレームに関連していない場合は、任意のRオブジェクトをファイル[^export]に保存することも可能です:
R
saveRDS(iris, file="/tmp/iris.rds")
ただし、このファイルの内容は人間が読める形式ではないため、Stack
Overflowに直接投稿することはできません。代わりに、readRDS()コマンドで閲覧可能な相手にメールで送信することができます(ここでは、ダウンロードしたファイルがユーザーのホームディレクトリ内のDownloadsフォルダにあることを前提としています):
R
some_data <- readRDS(file="~/Downloads/iris.rds")
最後になりますが、同様に重要な点として、必ずsessionInfo()の出力結果を記載してください。そこにはお使いのプラットフォームやRおよびパッケージのバージョン、その他問題の把握に非常に役立つ情報が含まれているためです。
R
sessionInfo()
出力
R version 4.6.2 (2026-10-29)
Platform: aarch64-apple-darwin23
Running under: macOS Golden Gate 27.2
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
time zone: Asia/Tokyo
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_4.6.2 cli_3.6.6 tools_4.6.2 otel_0.2.0 knitr_1.52
[6] xfun_0.62 rlang_1.3.0 png_0.1-9 evaluate_1.0.5
どこで助けを求めればよいでしょうか?
- コース中に隣に座っている方です。ワークショップ中は、遠慮なく隣の方と話をしたり、解答を照らし合わせたり、助けを求めたりしてください。
- 親切な仲間:自分より経験豊富な知り合いがいれば、快く助けてくれるかもしれません。
- Stack Overflow:質問にまだ回答がなく、かつ適切に構成されていれば、5分以内に回答が得られる可能性が高いです。良い質問の仕方に関するガイドラインに従うようにしてください。
- R-helpメーリングリスト:非常に多くの人が閲覧しており(Rコアチームの大半も含まれます)、投稿数も多いですが、口調がかなり堅苦しく、新規ユーザーを必ずしも温かく迎えてくれるとは限りません。質問内容が適切であれば、おそらく非常に早く回答が得られるでしょうが、絵文字付きの返信を期待しないでください。また、こここそ他のどこよりも、正しい用語を使うよう心がけてください(そうしないと、質問に答えるのではなく、言葉の誤用を指摘する回答が返ってくる可能性があります)。さらに、特定のパッケージではなくbase関数に関する質問の方が、よりうまく回答を得られるでしょう。
- 特定のパッケージに関する質問の場合は、そのパッケージのメーリングリストがないか確認してください。通常、
packageDescription("name-of-package")でアクセスできるパッケージの DESCRIPTION ファイルに記載されています。また、パッケージの作者に直接メールを送るか、コードリポジトリ(例:GitHub)で issue を起票してみるのもよいでしょう。 - また、特定のトピックに特化したメーリングリスト(GIS、系統学など)もあり、全リストはR-projectのメーリングリストページでご確認いただけます。
Rパッケージ
パッケージの読み込み
前述の通り、RではRパッケージが極めて重要な役割を果たします。パッケージがインストールされていることを前提として、その機能を利用するには、まずパッケージをロードする必要があります。これはlibrary()関数を用いて行います。以下では、ggplot2をロードします。
R
library("ggplot2")
パッケージのインストール
デフォルトのパッケージリポジトリは Comprehensive R Archive
Network (CRAN) であり、CRANで提供されているパッケージはすべて
install.packages()
関数でインストールできます。以下では、例えば後ほど学習する
dplyr パッケージをインストールします。
R
install.packages("dplyr")
このコマンドを実行すると、dplyrパッケージとそのすべての依存関係、すなわち動作に必要となるすべてのパッケージがインストールされます。
もう一つの主要なRパッケージリポジトリは、Bioconductorによって管理されています。Bioconductorパッケージは、専用のパッケージであるBiocManagerを用いて管理・インストールされます。このパッケージは、CRANから以下のようにインストールできます。
R
install.packages("BiocManager")
SummarizedExperiment(後ほど使用します)やDESeq2(RNA-Seq解析用)、その他BioconductorやCRANの各種パッケージは、BiocManager::installでインストールできます。
R
BiocManager::install("SummarizedExperiment")
BiocManager::install("DESeq2")
デフォルトでは、BiocManager::install()はインストール済みのすべてのパッケージを確認し、より新しいバージョンがないかもチェックします。新しいバージョンがある場合は、それを表示してUpdate all/some/none? [a/s/n]:するかどうかを尋ね、回答を待ちます。可能な限り最新のパッケージバージョンを維持することが望ましいですが、実際には、パッケージを読み込む前の新規Rセッションでのみパッケージを更新することを推奨します。
- RとRStudioを使い始めましょう
コマンドラインコンソールから直接Rを使用する場合とは異なります。Rと連携・統合する他のソフトウェアも存在しますが、RStudioは初心者に特に適している一方で、非常に高度な機能を数多く備えています。↩︎
すなわち、バイオインフォマティクスのデータ解析など、Rに新たな機能を付与するアドオンのことです。↩︎
本コースでは、バイオインフォマティクスを、生物学的または生物医学的データに適用されるデータサイエンスとして捉えます。↩︎
分析を実行してから数ヶ月、あるいは数年後のあなたの姿が、まさにその姿である可能性があり、また実際にそうなる可能性も非常に高いでしょう。↩︎
ここでは(統計学を除き)その大部分を紹介しますが、Rで実現できる膨大な可能性のほんの一端に触れるにとどまります。↩︎
Content from R入門
最終更新日:2026-11-07 | ページの編集
所要時間: 120分
概要
質問
- Rでの最初のコマンド
目的
- Rに関連する以下の用語を定義してください。オブジェクト、代入、呼び出し、関数、引数、オプション。
- Rでオブジェクトに値を代入します。
- オブジェクトの命名方法を学びましょう。
- コメントを用いてスクリプトに説明を添えてください。
- Rで簡単な算術演算を行います。
- 関数を呼び出し、引数を用いてデフォルトオプションを変更します。
- ベクトルの内容を確認し、その内容を操作します。
- ベクトルから部分集合を抽出し、値を取得します。
- 欠損値を含むベクトルを解析します。
このエピソードは、Data Carpentriesの「生態学者のためのRによるデータ分析と可視化」レッスンに基づいています。
Rでのオブジェクトの作成
コンソールに数式を入力するだけで、Rから出力を得ることができます。
R
3 + 5
出力
[1] 8
R
12 / 7
出力
[1] 1.714286
しかし、有用で興味深いことを行うには、オブジェクトに値を代入する必要があります。オブジェクトを作成するには、名前に続けて代入演算子
<- と、代入したい値を指定します:
R
weight_kg <- 55
<-
は代入演算子です。右側の値を左側のオブジェクトに代入します。したがって、x <- 3
を実行した後、x の値は 3
となります。この矢印は、3 が x
に入る
と読み解くことができます。歴史的な経緯から、代入に =
を使用することもできますが、あらゆる文脈で使用できるわけではありません。構文にわずかな違いがあるため、代入には常に
<- を使用するのが望ましい慣行です。
RStudioでは、Alt +
-(-キーと同時にAltを押す)と入力すると、PCでは1回のキー操作で<-が入力されます。一方、MacではOption
+
-(-キーと同時にOptionを押す)と入力することで同様の操作が可能です。
変数の命名
オブジェクトには、x、current_temperature、subject_id
などの任意の名前を付けられます。オブジェクト名は明示的で、長すぎないようにしてください。数字で始めることはできません(2x
は無効ですが、x2
は有効です)。Rは大文字と小文字を区別します(例:weight_kg
は Weight_kg
とは異なります)。Rの基本関数の名前であるため、使用できない名前もあります(例:if、else、for。完全なリストについては
Rの予約語
のマニュアルページを参照してください)。一般的に、許可されていても、他の関数名(例:c、T、mean、data、df、weights)は使用しないのが最善です。判断に迷う場合は、ヘルプを確認してその名前が既に使用されていないか調べてください。また、my.dataset
のようにオブジェクト名の中にドット(.)を含めることも避けるのが最善です。歴史的な理由から、Rには名前にドットを含む関数が多く存在しますが、ドットはR(メソッド用)や他のプログラミング言語において特別な意味を持つため、避けるのが最善です。また、オブジェクト名には名詞を、関数名には動詞を使用することが推奨されます。コードのスタイル(スペースの配置やオブジェクトの命名方法など)を一貫させることが重要です。一貫したコーディングスタイルを採用することで、将来の自分や共同作業者にとってコードがより読みやすくなります。Rでは、Google
のスタイルガイド、tidyverse
のスタイル、および Bioconductor
スタイルガイド などが一般的です。tidyverse
のガイドは非常に包括的で、最初は圧倒されるかもしれません。lintr
パッケージをインストールすれば、コードのスタイルに関する問題を自動的にチェックできます。
オブジェクトと変数:
Rでobjectsと呼ばれるものは、他の多くのプログラミング言語ではvariablesと呼ばれます。文脈によっては、objectとvariableは大きく異なる意味を持つことがあります。ただし、本レッスンでは、これら2つの用語は同義語として使用されます。詳細については、こちらをご覧ください。
オブジェクトに値を代入する際、Rは何も出力しません。括弧を使うか、オブジェクト名を入力することで、Rに値を強制的に出力させることができます。
R
weight_kg <- 55 # doesn't print anything
(weight_kg <- 55) # but putting parenthesis around the call prints the value of `weight_kg`
出力
[1] 55
R
weight_kg # and so does typing the name of the object
出力
[1] 55
Rのメモリにweight_kgが格納されたので、これを使って計算を行うことができます。例えば、この重量をポンドに変換したい場合が考えられます(ポンドでの重量はkgでの重量の2.2倍です):
R
2.2 * weight_kg
出力
[1] 121
また、オブジェクトに新しい値を代入することで、その値を変更することもできます。
R
weight_kg <- 57.5
2.2 * weight_kg
出力
[1] 126.5
つまり、あるオブジェクトに値を代入しても、他のオブジェクトの値は変わりません。例えば、動物の体重(ポンド)を新しいオブジェクト
weight_lb に格納してみましょう。
R
weight_lb <- 2.2 * weight_kg
次に、weight_kgを100に変更してください。
R
weight_kg <- 100
課題:
オブジェクト weight_lb
の現在の内容は、何だと思いますか? 126.5 でしょうか、それとも 220
でしょうか?
コメント
Rにおけるコメント記号は#です。スクリプト内で#より右側にある記述は、Rによって無視されます。スクリプトにメモや説明を残しておくと便利です。
RStudioでは、段落のコメントアウトやコメント解除を簡単に行えます。コメントアウトしたい行を選択した後、キーボードの Ctrl + Shift + C を同時に押してください。1行だけをコメントアウトしたい場合は、その行の任意の位置にカーソルを置き(つまり、行全体を選択する必要はありません)、Ctrl + Shift + C を押してください。
課題
以下において、各文の後の値は何になりますか?
R
mass <- 47.5 # mass?
age <- 122 # age?
mass <- mass * 2.0 # mass?
age <- age - 20 # age?
mass_index <- mass/age # mass_index?
関数とその引数
関数とは、代入操作などを含むより複雑な一連のコマンドを自動化する「定型スクリプト」のことです。多くの関数はあらかじめ定義されているか、あるいはRのパッケージをインポートすることで利用可能になります(これについては後述します)。関数は通常、引数と呼ばれる1つ以上の入力を受け取ります。関数は多くの場合(必ずしもそうとは限りませんが)、値を返します。典型的な例は関数
sqrt()
です。入力(引数)は数値である必要があり、戻り値(実質的には出力)はその数値の平方根となります。関数を実行すること(「動かす」こと)を、関数を呼び出すと言います。関数呼び出しの例は以下の通りです:
R
b <- sqrt(a)
ここでは、aの値がsqrt()関数に渡され、sqrt()関数が平方根を計算してその値を返し、それがオブジェクトbに代入されます。この関数は引数を1つしか取らないため、非常にシンプルです。
関数の戻り値「value」は、必ずしも数値(sqrt()のように)である必要はありません。また、単一の要素である必要もなく、複数の要素の集合や、データセットであることも可能です。これについては、データファイルをRに読み込む際に見ていきましょう。
引数には、数値やファイル名だけでなく、その他のオブジェクトも含め、あらゆるものを指定できます。各引数が具体的に何を意味するかは関数によって異なり、ドキュメント(以下参照)で確認する必要があります。一部の関数では、ユーザーが指定できる引数のほか、指定しなかった場合にデフォルト値が適用される引数があります。これらは「オプション」と呼ばれます。オプションは通常、「不正な値」を無視するかどうかや、プロットで使用する記号など、関数の動作を変更するために用いられます。ただし、特定の値を指定したい場合は、デフォルト値の代わりに任意の値を指定することができます。
複数の引数を受け取れる関数を試してみましょう。round()。
R
round(3.14159)
出力
[1] 3
ここでは、round()を引数3.14159のみで呼び出しており、3という値が返されました。これは、デフォルトで最近接の整数に丸められるためです。より多くの桁数が必要な場合は、round関数に関する情報を取得することでその方法を確認できます。args(round)を使用するか、?roundを使ってこの関数のヘルプを参照することができます。
R
args(round)
出力
function (x, digits = 0, ...)
NULL
R
?round
桁数を変更したい場合は、digits=2
や任意の桁数を入力すればよいことがわかります。
R
round(3.14159, digits = 2)
出力
[1] 3.14
引数を定義時と全く同じ順序で指定すれば、名前を付ける必要はありません。
R
round(3.14159, 2)
出力
[1] 3.14
また、引数に名前を付けておけば、その順序を入れ替えることもできます。
R
round(digits = 2, x = 3.14159)
出力
[1] 3.14
関数呼び出しの際には、必須の引数(丸め処理を行う数値など)を最初に配置し、すべてのオプション引数の名前を指定するのが良い習慣です。そうしないと、コードを読む人が、馴染みのない引数を持つ関数の定義を調べて、何を行っているのかを理解しなければならない場合があります。また、引数の名前を指定することで、将来的に関数インターフェースに変更が加えられ、既存の引数の間に新しい引数が追加される可能性に対する備えにもなります。
ベクトルとデータ型
ベクトルはRで最も一般的かつ基本的なデータ型であり、実質的にRの主力となるデータ型です。ベクトルは、数値や文字などの一連の値で構成されます。c()関数を使用して、一連の値をベクトルに代入できます。例えば、動物の体重のベクトルを作成し、それを新しいオブジェクトweight_gに代入することができます:
R
weight_g <- c(50, 60, 65, 82)
weight_g
出力
[1] 50 60 65 82
ベクトルには文字を含めることもできます。
R
molecules <- c("dna", "rna", "protein")
molecules
出力
[1] "dna" "rna" "protein"
ここでは「dna」や「rna」などの引用符が不可欠です。引用符がないと、Rはdna、rna、proteinというオブジェクトが存在すると判断します。これらのオブジェクトはRのメモリ上に存在しないため、エラーメッセージが表示されます。
ベクトルの内容を確認できる関数が数多く用意されています。length()
を使うと、特定のベクトルにいくつの要素が含まれているかが分かります。
R
length(weight_g)
出力
[1] 4
R
length(molecules)
出力
[1] 3
ベクトルの重要な特徴は、すべての要素が同じデータ型であることです。関数
class()
は、オブジェクトのクラス(要素の型)を示します。
R
class(weight_g)
出力
[1] "numeric"
R
class(molecules)
出力
[1] "character"
str()関数は、オブジェクトの構造とその要素の概要を提供します。大規模で複雑なオブジェクトを扱う際に有用な関数です。
R
str(weight_g)
出力
num [1:4] 50 60 65 82
R
str(molecules)
出力
chr [1:3] "dna" "rna" "protein"
c()関数を使えば、ベクトルに他の要素を追加できます。
R
weight_g <- c(weight_g, 90) # add to the end of the vector
weight_g <- c(30, weight_g) # add to the beginning of the vector
weight_g
出力
[1] 30 50 60 65 82 90
1行目では、元のベクトル weight_g の末尾に値
90 を追加し、その結果を weight_g
に保存します。次に、先頭に値 30 を追加し、同様にその結果を
weight_g に保存します。
これを繰り返すことで、ベクトルを拡張したり、データセットを構築したりできます。プログラミングを進める中で、収集や計算を行っている結果を追加する際にも役立つでしょう。
アトミックベクトルはRで最も単純なデータ型であり、単一の型からなる線形ベクトルです。上記では、Rが使用する主要な6種類のアトミックベクトル型のうち、"character"と"numeric"(または"double")の2つを確認しました。これらは、すべてのRオブジェクトを構成する基本的な構成要素です。その他の4種類のアトミックベクトル型は以下の通りです。
-
TRUEおよびFALSE用の"logical"(論理型データ) - 整数には
"integer"を使用します(例:2L。ここでLはRに対してそれが整数であることを示します)。 -
"complex"は実部と虚部を持つ複素数(例:1 + 4i)を表します。これらについては以上です。 -
"raw"以降では触れないビットストリーム用
typeof()関数にベクトルを引数として渡すことで、その型を確認できます。
ベクトルは、Rが使用する数多くのデータ構造のひとつです。その他の重要なものには、リスト
(list)、行列 (matrix)、データフレーム
(data.frame)、因子 (factor)、配列
(array) があります。
課題:
アトミックベクトルには、文字型、数値型(またはdouble型)、整数型、論理型があることを確認しました。では、これらの異なる型を一つのベクトルに混在させようとすると、どうなるでしょうか。
Rはこれらを暗黙的にすべて同じ型に変換します。
課題:
それぞれの例で何が起こるでしょうか?(ヒント:class()
を使ってオブジェクトのデータ型を確認し、その名前を入力して何が起こるか確かめてみてください):
R
num_char <- c(1, 2, 3, "a")
num_logical <- c(1, 2, 3, TRUE, FALSE)
char_logical <- c("a", "b", "c", TRUE)
tricky <- c(1, 2, 3, "4")
R
class(num_char)
出力
[1] "character"
R
num_char
出力
[1] "1" "2" "3" "a"
R
class(num_logical)
出力
[1] "numeric"
R
num_logical
出力
[1] 1 2 3 1 0
R
class(char_logical)
出力
[1] "character"
R
char_logical
出力
[1] "a" "b" "c" "TRUE"
R
class(tricky)
出力
[1] "character"
R
tricky
出力
[1] "1" "2" "3" "4"
課題:
なぜそうなると思いますか?
ベクトルは単一のデータ型しか持つことができません。Rは、このベクトルの内容を変換(coerce)し、情報を一切損なわない共通の型を見つけ出そうとします。
課題:
次の例において、combined_logical
に含まれる値のうち、"TRUE"(文字列として)となるものはいくつありますか?
R
num_logical <- c(1, 2, 3, TRUE)
char_logical <- c("a", "b", "c", TRUE)
combined_logical <- c(num_logical, char_logical)
1つだけです。過去のデータ型の記憶はなく、型変換はベクトルが最初に評価される際に行われます。したがって、num_logical
内の TRUE は、combined_logical で
"1" に変換される前に、まず 1
に変換されます。
R
combined_logical
出力
[1] "1" "2" "3" "1" "a" "b" "c" "TRUE"
課題:
Rでは、あるクラスのオブジェクトを別のクラスに変換することを「強制変換」と呼びます。これらの変換は階層構造に基づいて行われ、一部の型は他の型へ優先的に強制変換されます。これらのデータ型がどのように強制変換されるか、その階層構造を示す図を描いていただけますか?
論理型 → 数値型 → 文字列型 ← 論理型
ベクトルの部分抽出
ベクトルから1つまたは複数の値を抽出したい場合は、角括弧で1つまたは複数のインデックスを指定する必要があります。例えば:
R
molecules <- c("dna", "rna", "peptide", "protein")
molecules[2]
出力
[1] "rna"
R
molecules[c(3, 2)]
出力
[1] "peptide" "rna"
また、インデックスを繰り返すことで、元のオブジェクトよりも要素数の多いオブジェクトを作成することもできます。
R
more_molecules <- molecules[c(1, 2, 3, 2, 1, 4)]
more_molecules
出力
[1] "dna" "rna" "peptide" "rna" "dna" "protein"
Rのインデックスは1から始まります。Fortran、MATLAB、Julia、Rなどのプログラミング言語が1から数えるのは、それが人間にとって一般的だからです。C系言語(C++、Java、Perl、Pythonを含む)が0から数えるのは、コンピュータにとってその方が処理が簡単だからです。
最後に、負のインデックスを用いることで、ベクトルの指定した要素を除くすべての要素を取得することも可能です。
R
molecules ## all molecules
出力
[1] "dna" "rna" "peptide" "protein"
R
molecules[-1] ## all but the first one
出力
[1] "rna" "peptide" "protein"
R
molecules[-c(1, 3)] ## all but 1st/3rd ones
出力
[1] "rna" "protein"
R
molecules[c(-1, -3)] ## all but 1st/3rd ones
出力
[1] "rna" "protein"
条件付きサブセット抽出
サブセットを抽出するもう一つの一般的な方法は、論理ベクトルを使用することです。TRUEは同じインデックスを持つ要素を選択しますが、FALSEは選択しません。
R
weight_g <- c(21, 34, 39, 54, 55)
weight_g[c(TRUE, FALSE, TRUE, TRUE, FALSE)]
出力
[1] 21 39 54
通常、これらの論理ベクトルは手動で入力されるものではなく、他の関数や論理判定の結果として得られるものです。例えば、50を超える値のみを選択したい場合は、次のようにします。
R
## will return logicals with TRUE for the indices that meet
## the condition
weight_g > 50
出力
[1] FALSE FALSE FALSE TRUE TRUE
R
## so we can use this to select only the values above 50
weight_g[weight_g > 50]
出力
[1] 54 55
&(両方の条件が真、AND)または|(少なくとも一方の条件が真、OR)を使用して、複数のテストを組み合わせることができます。
R
weight_g[weight_g < 30 | weight_g > 50]
出力
[1] 21 54 55
R
weight_g[weight_g >= 30 & weight_g == 21]
出力
numeric(0)
ここで、<は「未満」、>は「より大きい」、>=は「以上」、==は「等しい」を表します。二重等号==は左辺と右辺の数値的な等価性を判定するもので、変数への代入を行う単一の=記号(<-と同様)と混同しないようご注意ください。
ベクトル内で特定の文字列を検索することは、よくある作業です。複数の値との等価性を判定するには「or」演算子
| を使用できますが、すぐに煩雑になります。関数
%in%
を使用すると、検索ベクトルの要素のいずれかが見つかるかどうかを判定できます:
R
molecules <- c("dna", "rna", "protein", "peptide")
molecules[molecules == "rna" | molecules == "dna"] # returns both rna and dna
出力
[1] "dna" "rna"
R
molecules %in% c("rna", "dna", "metabolite", "peptide", "glycerol")
出力
[1] TRUE TRUE FALSE TRUE
R
molecules[molecules %in% c("rna", "dna", "metabolite", "peptide", "glycerol")]
出力
[1] "dna" "rna" "peptide"
課題:
なぜ "four" > "five" が TRUE
を返すのか、分かりますか?
R
"four" > "five"
出力
[1] TRUE
文字列に対して > や <
を使用する場合、R はそれらのアルファベット順を比較します。ここでは
"four" が "five"
の後に来るため、前者の方が後者よりも大きいと判断されます。
名称
ベクトルの各要素に名前を付けることができます。以下のコードチャンクでは、名前が付与されていない初期のベクトルと、名前の設定方法および取得方法を示しています。
R
x <- c(1, 5, 3, 5, 10)
names(x) ## no names
出力
NULL
R
names(x) <- c("A", "B", "C", "D", "E")
names(x) ## now we have names
出力
[1] "A" "B" "C" "D" "E"
ベクトルに名前が付いている場合、インデックスに加えて名前で要素にアクセスすることができます。
R
x[c(1, 3)]
出力
A C
1 3
R
x[c("A", "C")]
出力
A C
1 3
欠損データ
Rはデータセットの解析を目的として設計されているため、欠損データという概念を備えています(他のプログラミング言語では一般的ではありません)。欠損データは、ベクトル内でNAとして表現されます。
数値演算を行う際、扱うデータに欠損値が含まれていると、ほとんどの関数はNAを返します。この仕様により、欠損値を扱うケースを見落とすことが難しくなります。引数na.rm = TRUEを指定することで、欠損値を無視して結果を算出できます。
R
heights <- c(2, 4, 4, NA, 6)
mean(heights)
出力
[1] NA
R
max(heights)
出力
[1] NA
R
mean(heights, na.rm = TRUE)
出力
[1] 4
R
max(heights, na.rm = TRUE)
出力
[1] 6
データに欠損値が含まれている場合は、is.na()、na.omit()、complete.cases()
の各関数に慣れておくとよいでしょう。例については以下をご覧ください。
R
## Extract those elements which are not missing values.
heights[!is.na(heights)]
出力
[1] 2 4 4 6
R
## Returns the object with incomplete cases removed.
## The returned object is an atomic vector of type `"numeric"`
## (or `"double"`).
na.omit(heights)
出力
[1] 2 4 4 6
attr(,"na.action")
[1] 4
attr(,"class")
[1] "omit"
R
## Extract those elements which are complete cases.
## The returned object is an atomic vector of type `"numeric"`
## (or `"double"`).
heights[complete.cases(heights)]
出力
[1] 2 4 4 6
課題:
- このインチ単位の身長のベクトルを用いて、NAを除去した新しいベクトルを作成してください。
R
heights <- c(63, 69, 60, 65, NA, 68, 61, 70, 61, 59, 64, 69, 63, 63, NA, 72, 65, 64, 70, 63, 65)
- 関数
median()を使用して、ベクトルheightsの中央値を計算してください。 - Rを使って、集合の中で67インチより背が高い人が何人いるか調べてください。
R
heights_no_na <- heights[!is.na(heights)]
## or
heights_no_na <- na.omit(heights)
R
median(heights, na.rm = TRUE)
出力
[1] 64
R
heights_above_67 <- heights_no_na[heights_no_na > 67]
length(heights_above_67)
出力
[1] 6
ベクトルの生成
コンストラクタ
さまざまな型のベクトルを生成するための関数がいくつか用意されています。数値型のベクトルを生成するには、出力ベクトルの長さを引数として指定するnumeric()コンストラクタを使用できます。値は0で初期化されます。
R
numeric(3)
出力
[1] 0 0 0
R
numeric(10)
出力
[1] 0 0 0 0 0 0 0 0 0 0
なお、長さ0の数値ベクトルを要求すると、まさにそれが得られます。
R
numeric(0)
出力
numeric(0)
文字型と論理型に対しても、それぞれcharacter()およびlogical()という同様のコンストラクタが用意されています。
課題:
文字型および論理型ベクトルのデフォルト値は何ですか?
R
character(2) ## the empty character
出力
[1] "" ""
R
logical(2) ## FALSE
出力
[1] FALSE FALSE
要素を複製します
rep
関数を使用すると、値を指定した回数だけ繰り返すことができます。例えば、値
-1 で長さ 5
の数値ベクトルを初期化したい場合は、次のように記述できます。
R
rep(-1, 5)
出力
[1] -1 -1 -1 -1 -1
同様に、収集するデータに前提条件を設けず、まずは欠損値を含むベクトルを生成する方法もよく用いられます。
R
rep(NA, 5)
出力
[1] NA NA NA NA NA
repは、任意の長さのベクトル(上記では長さ1のベクトルを使用しました)および任意の型を入力として受け付けます。例えば、値1、2、3を5回繰り返したい場合は、次のように記述します。
R
rep(c(1, 2, 3), 5)
出力
[1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
課題:
値 1、2、3 をそれぞれ 5 回繰り返しつつ、その順に 1 が 5 個、2 が 5
個、3 が 5 個得たい場合はどうすればよいでしょうか。方法は 2
通りあります。詳細は ?rep または ?sort
を参照してください。
R
rep(c(1, 2, 3), each = 5)
出力
[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3
R
sort(rep(c(1, 2, 3), 5))
出力
[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3
数列の生成
もう一つ非常に便利な関数に、数値列を生成するseqがあります。例えば、1から20まで2刻みで整数列を生成するには、次のように記述します。
R
seq(from = 1, to = 20, by = 2)
出力
[1] 1 3 5 7 9 11 13 15 17 19
byのデフォルト値は1です。1刻みで値を順次生成する操作は頻繁に用いられるため、ショートカットが用意されています。
R
seq(1, 5, 1)
出力
[1] 1 2 3 4 5
R
seq(1, 5) ## default by
出力
[1] 1 2 3 4 5
R
1:5
出力
[1] 1 2 3 4 5
1から20までの数値を、最終的に長さ3となるように生成するには、次のようにします。
R
seq(from = 1, to = 20, length.out = 3)
出力
[1] 1.0 10.5 20.0
無作為抽出と置換
最後に、有用な関数群として、ランダムデータを生成するものがあります。最初のsampleは、別のベクトルのランダムな順列を生成します。例えば、10人の学生の口頭試験のランダムな順序を決めるには、まず各学生に1から10までの番号を割り当て(例えば名前のアルファベット順に基づいて)、次に次のようにします:
R
sample(1:10)
出力
[1] 9 4 7 1 2 5 3 10 6 8
追加の引数を指定しない場合、sampleはベクトルの全要素の順列を返します。特定のサイズのランダムサンプルを取得したい場合は、この値を第2引数として指定します。以下では、あらかじめ定義されたlettersベクトルに含まれるアルファベットから5文字をランダムに抽出します。
R
sample(letters, 5)
出力
[1] "s" "a" "u" "x" "j"
入力ベクトルよりも大きな出力を得たい場合や、一部の要素を複数回抽出したい場合は、replace引数をTRUEに設定する必要があります:
R
sample(1:5, 10, replace = TRUE)
出力
[1] 2 1 5 5 1 1 5 5 2 2
課題:
上記の関数を試してみると、サンプルが確かにランダムに抽出されており、同じ順列が二度現れることはないことに気づいたはずです。こうしたランダム抽出を再現するには、ランダムサンプルを抽出する前に、set.seed()
を用いて乱数生成のシード値を手動で設定することができます。
隣の人と一緒にこの機能を試してみましょう。まず、1:10のランダムな順列を2つ独立して作成し、結果が異なることを確認してください。
次に、例えば set.seed(123)
でシード値を設定し、ランダム抽出を再度行います。これで、同じランダム抽出結果が得られることが確認できます。
別のシード値を設定して繰り返してください。
さまざまな順列
R
sample(1:10)
出力
[1] 9 1 4 3 6 2 5 8 10 7
R
sample(1:10)
出力
[1] 4 9 7 6 1 10 8 3 2 5
シード値 123 での同一の順列
R
set.seed(123)
sample(1:10)
出力
[1] 3 10 2 8 6 9 1 7 5 4
R
set.seed(123)
sample(1:10)
出力
[1] 3 10 2 8 6 9 1 7 5 4
別のシード値
R
set.seed(1)
sample(1:10)
出力
[1] 9 4 7 1 2 5 3 10 6 8
R
set.seed(1)
sample(1:10)
出力
[1] 9 4 7 1 2 5 3 10 6 8
正規分布からの標本抽出
最後に見ていく関数はrnormで、正規分布から無作為標本を抽出します。平均0と100、標準偏差1と5の2つの正規分布(N(0,
1)およびN(100, 5)と表記)を以下に示します。

3つの引数 n、mean、sd
は、標本サイズと、正規分布のパラメータ(すなわち平均とその標準偏差)を指定します。後者のデフォルト値は
0 と 1 です。
R
rnorm(5)
出力
[1] 0.69641761 0.05351568 -1.31028350 -2.12306606 -0.20807859
R
rnorm(5, 2, 2)
出力
[1] 1.3744268 -0.1164714 2.8344472 1.3690969 3.6510983
R
rnorm(5, 100, 5)
出力
[1] 106.45636 96.87448 95.62427 100.71678 107.12595
スクリプトの書き方とRのデータ構造の基礎を学んだところで、いよいよより大規模なデータを扱い、データフレームについて学んでいきましょう。
- Rの操作方法
Content from データから始める
最終更新日:2026-11-07 | ページの編集
所要時間: 60分
概要
質問
- Rによる初めてのデータ分析
目的
-
data.frameとは何か説明してください。 - .csvファイルから外部データをデータフレームに読み込みます。
- データフレームの内容を要約します。
- 因子とは何か説明してください。
- 文字列と因子を相互に変換します。
- 因子の順序を入れ替え、名称を変更します。
- 日付をフォーマットします。
- データをエクスポートして保存します。
このエピソードは、Data Carpentriesの「生態学者のためのRによるデータ分析と可視化」レッスンに基づいています。
遺伝子発現データの提示
Blackmore et al. (2017)が発表したデータの一部である『The effect of upper-respiratory infection on transcriptomic changes in the CNS』を使用します。本研究の目的は、上気道感染症が感染後の小脳および脊髄で起こるRNA転写の変化に及ぼす影響を明らかにすることでした。性別を一致させた8週齢のC57BL/6マウスに、生理食塩水またはインフルエンザA型を経鼻投与し、0日目(非感染)、4日目、および8日目にRNA-seqを用いて小脳および脊髄組織におけるトランスクリプトームの変化を評価しました。
データセットはカンマ区切り値(CSV)ファイルとして保存されています。各行には単一のRNA発現測定に関する情報が含まれており、最初の11列はそれぞれ以下を表しています。
| 列 | 説明 |
|---|---|
| gene | 測定された遺伝子名 |
| sample | 遺伝子発現を測定したサンプル名 |
| expression | 遺伝子発現の値 |
| organism | 生物体/種 - ここではすべてのデータがマウス由来です。 |
| age | マウスの年齢(ここではすべてのマウスが8週齢でした) |
| sex | マウスの性別 |
| infection | マウスの感染状態、すなわちインフルエンザA型への感染の有無です。 |
| strain | A型インフルエンザ株。 |
| time | 感染期間(日数)。 |
| tissue | 遺伝子発現実験に用いられた組織、すなわち小脳または脊髄。 |
| mouse | マウスの一意の識別子。 |
Rの関数download.file()を使用して遺伝子発現データを含むCSVファイルをダウンロードし、read.csv()を用いてCSVファイルの内容をクラスdata.frameのオブジェクトとしてメモリに読み込みます。download.fileコマンドの最初の項目は、ソースURLを示す文字列です。このソースURLは、GitHubリポジトリからCSVファイルをダウンロードします。カンマ("data/rnaseq.csv")以降のテキストは、ローカルマシン上のファイルの保存先です。ファイルをダウンロードするための"data"というフォルダがマシン上に必要になります。したがって、このコマンドはリモートファイルをダウンロードし、"rnaseq.csv"という名前を付けて、既存の"data"というフォルダに追加します。
R
download.file(url = "https://github.com/carpentries-incubator/bioc-intro/raw/main/episodes/data/rnaseq.csv",
destfile = "data/rnaseq.csv")
これでデータを読み込む準備が整いました。
R
rna <- read.csv("data/rnaseq.csv")
この記述では何も出力されません。ご記憶にある通り、代入を行っても何も表示されないためです。データが読み込まれたか確認したい場合は、データフレーム名を入力することでその内容を表示できます。
R
rna
うわあ……かなりの出力量ですね。少なくとも、データが正常に読み込まれたということですね。関数
head()
を使って、このデータフレームの先頭(最初の6行)を確認してみましょう:
R
head(rna)
出力
gene sample expression organism age sex infection strain time
1 Asl GSM2545336 1170 Mus musculus 8 Female InfluenzaA C57BL/6 8
2 Apod GSM2545336 36194 Mus musculus 8 Female InfluenzaA C57BL/6 8
3 Cyp2d22 GSM2545336 4060 Mus musculus 8 Female InfluenzaA C57BL/6 8
4 Klk6 GSM2545336 287 Mus musculus 8 Female InfluenzaA C57BL/6 8
5 Fcrls GSM2545336 85 Mus musculus 8 Female InfluenzaA C57BL/6 8
6 Slc2a4 GSM2545336 782 Mus musculus 8 Female InfluenzaA C57BL/6 8
tissue mouse ENTREZID
1 Cerebellum 14 109900
2 Cerebellum 14 11815
3 Cerebellum 14 56448
4 Cerebellum 14 19144
5 Cerebellum 14 80891
6 Cerebellum 14 20528
product
1 argininosuccinate lyase, transcript variant X1
2 apolipoprotein D, transcript variant 3
3 cytochrome P450, family 2, subfamily d, polypeptide 22, transcript variant 2
4 kallikrein related-peptidase 6, transcript variant 2
5 Fc receptor-like S, scavenger receptor, transcript variant X1
6 solute carrier family 2 (facilitated glucose transporter), member 4
ensembl_gene_id external_synonym chromosome_name gene_biotype
1 ENSMUSG00000025533 2510006M18Rik 5 protein_coding
2 ENSMUSG00000022548 <NA> 16 protein_coding
3 ENSMUSG00000061740 2D22 15 protein_coding
4 ENSMUSG00000050063 Bssp 7 protein_coding
5 ENSMUSG00000015852 2810439C17Rik 3 protein_coding
6 ENSMUSG00000018566 Glut-4 11 protein_coding
phenotype_description
1 abnormal circulating amino acid level
2 abnormal lipid homeostasis
3 abnormal skin morphology
4 abnormal cytokine level
5 decreased CD8-positive alpha-beta T cell number
6 abnormal circulating glucose level
hsapiens_homolog_associated_gene_name
1 ASL
2 APOD
3 CYP2D6
4 KLK6
5 FCRL2
6 SLC2A4
R
## Try also
## View(rna)
注
read.csv()
はフィールドがカンマで区切られていることを前提としていますが、一部の国ではカンマが小数点の区切り文字として、セミコロン
(;) がフィールドの区切り文字として使用されています。R
でこのような形式のファイルを読み込む場合は、read.csv2()
関数を使用できます。この関数は read.csv()
と全く同様に動作しますが、小数点およびフィールドの区切り文字に異なるパラメータを使用します。別の形式を扱う場合は、どちらもユーザーが指定可能です。詳細については、?read.csv
と入力して read.csv()
のヘルプを参照してください。また、タブ区切りのデータファイルを読み込むための
read.delim()
関数もあります。重要な点として、これらの関数はすべて、実際には異なる引数を持つメインの
read.table()
関数のラッパー関数となっています。そのため、上記のデータは
, のように区切り文字引数を指定して
read.table()
を使用して読み込むこともできました。コードは以下の通りです。
R
rna <- read.table(file = "data/rnaseq.csv",
sep = ",",
header = TRUE)
デフォルトでは read.table() の header 引数が FALSE
に設定されているため、ヘッダーを読み込むには header 引数を TRUE
に設定する必要があります。
データフレームとは何ですか?
データフレームは、ほとんどの表形式データにおける事実上のデータ構造であり、統計処理やプロットに用いられます。
データフレームは手動で作成することもできますが、一般的にはread.csv()やread.table()といった関数を用いて生成されます。言い換えれば、ハードディスク(またはウェブ)からスプレッドシートをインポートする際に作成されます。
データフレームとは、列がすべて同じ長さのベクトルである表形式のデータ表現のことです。列はベクトルであるため、各列には単一のデータ型(例:文字、整数、因子)が含まれている必要があります。例えば、以下は数値ベクトル、文字ベクトル、および論理ベクトルからなるデータフレームを示した図です。
これは、str()関数を用いてデータフレームのstr構造を確認することで分かります。
R
str(rna)
出力
'data.frame': 32428 obs. of 19 variables:
$ gene : chr "Asl" "Apod" "Cyp2d22" "Klk6" ...
$ sample : chr "GSM2545336" "GSM2545336" "GSM2545336" "GSM2545336" ...
$ expression : int 1170 36194 4060 287 85 782 1619 288 43217 1071 ...
$ organism : chr "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" ...
$ age : int 8 8 8 8 8 8 8 8 8 8 ...
$ sex : chr "Female" "Female" "Female" "Female" ...
$ infection : chr "InfluenzaA" "InfluenzaA" "InfluenzaA" "InfluenzaA" ...
$ strain : chr "C57BL/6" "C57BL/6" "C57BL/6" "C57BL/6" ...
$ time : int 8 8 8 8 8 8 8 8 8 8 ...
$ tissue : chr "Cerebellum" "Cerebellum" "Cerebellum" "Cerebellum" ...
$ mouse : int 14 14 14 14 14 14 14 14 14 14 ...
$ ENTREZID : int 109900 11815 56448 19144 80891 20528 97827 118454 18823 14696 ...
$ product : chr "argininosuccinate lyase, transcript variant X1" "apolipoprotein D, transcript variant 3" "cytochrome P450, family 2, subfamily d, polypeptide 22, transcript variant 2" "kallikrein related-peptidase 6, transcript variant 2" ...
$ ensembl_gene_id : chr "ENSMUSG00000025533" "ENSMUSG00000022548" "ENSMUSG00000061740" "ENSMUSG00000050063" ...
$ external_synonym : chr "2510006M18Rik" NA "2D22" "Bssp" ...
$ chromosome_name : chr "5" "16" "15" "7" ...
$ gene_biotype : chr "protein_coding" "protein_coding" "protein_coding" "protein_coding" ...
$ phenotype_description : chr "abnormal circulating amino acid level" "abnormal lipid homeostasis" "abnormal skin morphology" "abnormal cytokine level" ...
$ hsapiens_homolog_associated_gene_name: chr "ASL" "APOD" "CYP2D6" "KLK6" ...
data.frame オブジェクトの確認
データフレームの内容や構造を確認する際に、head() や
str()
といった関数がどのように役立つかはすでに見てきました。ここでは、データの内容や構造を把握するための関数を、網羅的ではありませんがいくつか紹介します。さっそく試してみましょう!
サイズ:
-
dim(rna)- 行数を第1要素、列数を第2要素とするベクトル(オブジェクトの次元)を返します。 -
nrow(rna)- 行数を返します。 -
ncol(rna)- 列数を返します。
内容:
-
head(rna)- 最初の6行を表示します。 -
tail(rna)- 直近の6行を表示します。
名称:
-
names(rna)- 列名を返します(data.frameオブジェクトの場合はcolnames()と同義です)。 -
rownames(rna)- 行名を返します。
概要:
-
str(rna)- オブジェクトの構造、およびクラス、各列の長さと内容に関する情報。 -
summary(rna)- 各列の要約統計量。
注:これらの関数の多くは「汎用的」であり、data.frame以外の種類のオブジェクトにも使用できます。
課題:
str(rna)の出力結果をもとに、以下の質問に答えられますか?
- オブジェクト
rnaのクラスは何ですか? - このオブジェクトの行数と列数はそれぞれいくつですか?
- クラス: データフレーム
- 行数: 32428、列数: 19
データフレームのインデックス指定と部分抽出
今回のrnaデータフレームには行と列(2次元)があります。ここから特定のデータを抽出したい場合は、目的の「座標」を指定する必要があります。まず行番号、次に列番号の順で指定します。ただし、これらの座標の指定方法によって、得られる結果のクラスが異なる点に注意してください。
R
# first element in the first column of the data frame (as a vector)
rna[1, 1]
# first element in the 6th column (as a vector)
rna[1, 6]
# first column of the data frame (as a vector)
rna[, 1]
# first column of the data frame (as a data.frame)
rna[1]
# first three elements in the 7th column (as a vector)
rna[1:3, 7]
# the 3rd row of the data frame (as a data.frame)
rna[3, ]
# equivalent to head_rna <- head(rna)
head_rna <- rna[1:6, ]
head_rna
:は、整数の数値ベクトルを昇順または降順に生成する特殊な関数です。例として1:10や10:1をお試しください。詳細は「ベクトルの生成」セクションを参照してください。
また、「-」記号を使用して、データフレームの特定のインデックスを除外することもできます。
R
rna[, -1] ## The whole data frame, except the first column
rna[-c(7:66465), ] ## Equivalent to head(rna)
データフレームは、インデックスを指定する(前述の通り)ほか、列名を直接指定することでもサブセットを抽出できます。
R
rna["gene"] # Result is a data.frame
rna[, "gene"] # Result is a vector
rna[["gene"]] # Result is a vector
rna$gene # Result is a vector
RStudioでは、自動補完機能を使って、列の正式名称を正しく取得できます。
課題
rnaデータセットの200行目のデータのみを含むdata.frame(rna_200)を作成してください。nrow()がdata.frameの行数を返したことに注目してください。
その数値を使って、初期の
rnaデータフレームから最後の行だけを抽出してください。それを
tail()を用いて最終行として表示されている内容と比較し、想定通りになっているか確認してください。行番号の代わりに
nrow()を使用して、最後の行を抽出してください。その最終行から新しいデータフレーム(
rna_last)を作成します。
nrow()を使用して、rnaデータフレームの中央にある行を抽出してください。この行の内容をrna_middleという名前のオブジェクトに格納してください。nrow()を上記の-表記と組み合わせ、rnaデータセットの1行目から6行目のみを保持した状態でhead(rna)の動作を再現してください。
R
## 1.
rna_200 <- rna[200, ]
## 2.
## Saving `n_rows` to improve readability and reduce duplication
n_rows <- nrow(rna)
rna_last <- rna[n_rows, ]
## 3.
rna_middle <- rna[n_rows / 2, ]
## 4.
rna_head <- rna[-(7:n_rows), ]
因子
因子はカテゴリカルデータを表します。これらはラベルに紐付けられた整数として格納され、順序付きまたは順序なしのいずれかとなります。因子は見た目(および多くの場合の挙動)は文字ベクトルのようですが、実際にはRでは整数ベクトルとして扱われます。そのため、文字列として扱う際には十分な注意が必要です。
一度作成されると、因子には水準と呼ばれるあらかじめ定義された値の集合しか含めることができません。デフォルトでは、Rは常に水準をアルファベット順に並べ替えます。例えば、2つの水準を持つ因子がある場合:
R
sex <- factor(c("male", "female", "female", "male", "female"))
Rは、レベル"female"に1を、レベル"male"に2を割り当てます(このベクトルの最初の要素は"male"ですが、fがmより前にあるためです)。これは関数levels()を使用することで確認でき、レベル数はnlevels()で調べることができます:
R
levels(sex)
出力
[1] "female" "male"
R
nlevels(sex)
出力
[1] 2
因子の順序が重要でない場合もあれば、意味があるため(例:「低」、「中」、「高」)、可視化を向上させるため、あるいは特定の種類の解析で必要であるため、順序を指定したい場合もあります。ここでは、sex
ベクトル内のレベルを並べ替える一つの方法を以下に示します。
R
sex ## current order
出力
[1] male female female male female
Levels: female male
R
sex <- factor(sex, levels = c("male", "female"))
sex ## after re-ordering
出力
[1] male female female male female
Levels: male female
Rのメモリ上では、これらの因子は整数(1, 2,
3)で表現されますが、因子は自己記述的であるため、整数よりも多くの情報を含んでいます。"female"や"male"は、1や2よりも記述的です。どれが「オス」なのでしょうか。整数データだけでは判断できません。一方、因子はこの情報が組み込まれています。これは、レベル数が多い場合(本例のデータセットにおける遺伝子バイオタイプなど)に特に有用です。
データが因子として格納されている場合、plot()関数を使えば、各因子水準に含まれる観測数を一目で確認できます。それでは、データ内のオスとメスの数を見てみましょう。
R
plot(sex)

文字列への変換
因子を文字ベクトルに変換するには、as.character(x)を使用します。
R
as.character(sex)
出力
[1] "male" "female" "female" "male" "female"
因子の名称変更
これらの因子の名前を変更したい場合は、その水準を変更するだけで十分です。
R
levels(sex)
出力
[1] "male" "female"
R
levels(sex) <- c("M", "F")
sex
出力
[1] M F F M F
Levels: M F
R
plot(sex)

課題:
- 「F」と「M」をそれぞれ「Female」と「Male」にリネームしてください。
R
levels(sex)
出力
[1] "M" "F"
R
levels(sex) <- c("Male", "Female")
- 動物名の引用符の欠落
- 「feel」列に1件の欠損があります(おそらく毛のある動物のいずれかのものと思われます)。
- weight列のカンマが1つ抜けています
課題:
以下の例の各列のクラスを予測できますか?
str(country_climate)を使って予想を確認してください。
期待通りでしたか? なぜですか? なぜそうではないのでしょうか?
データフレームを作成する際、最後の変数の後に
stringsAsFactors = TRUEを追加して再度試してみてください。今度はどうなりますか?read.csv()を使用してテキスト形式のスプレッドシートを R に読み込む際にも、stringsAsFactorsを設定できます。
R
country_climate <- data.frame(
country = c("Canada", "Panama", "South Africa", "Australia"),
climate = c("cold", "hot", "temperate", "hot/temperate"),
temperature = c(10, 30, 18, "15"),
northern_hemisphere = c(TRUE, TRUE, FALSE, "FALSE"),
has_kangaroo = c(FALSE, FALSE, FALSE, 1)
)
R
country_climate <- data.frame(
country = c("Canada", "Panama", "South Africa", "Australia"),
climate = c("cold", "hot", "temperate", "hot/temperate"),
temperature = c(10, 30, 18, "15"),
northern_hemisphere = c(TRUE, TRUE, FALSE, "FALSE"),
has_kangaroo = c(FALSE, FALSE, FALSE, 1)
)
str(country_climate)
出力
'data.frame': 4 obs. of 5 variables:
$ country : chr "Canada" "Panama" "South Africa" "Australia"
$ climate : chr "cold" "hot" "temperate" "hot/temperate"
$ temperature : chr "10" "30" "18" "15"
$ northern_hemisphere: chr "TRUE" "TRUE" "FALSE" "FALSE"
$ has_kangaroo : num 0 0 0 1
データ型の自動変換は、便利な場合もあれば、厄介な場合もあります。この機能が存在することを認識し、ルールを把握した上で、Rにインポートしたデータがデータフレーム内で正しい型になっているかを再確認してください。もしそうでない場合は、それを活用して、データ入力時に混入した可能性のある誤り(例えば、数値のみであるべき列に文字が入っているなど)を検出してください。
詳細については、こちらのRStudioチュートリアルをご覧ください。
行列
先に進む前に、データフレームについて学んだところで、パッケージのインストールについて復習し、新しいデータ型であるmatrixについて学びましょう。data.frameと同様に、行列には行と列という2つの次元があります。しかし、大きな違いは、matrixのすべてのセルがnumeric、character、logical…といった同一の型でなければならない点です。その点において、行列はdata.frameよりもvectorに近いと言えます。
行列のデフォルトコンストラクタはmatrixです。行列に格納する値のベクトルと、行数および/または列数1を指定します。値は、以下に示すように列方向に配置されます。
R
m <- matrix(1:9, ncol = 3, nrow = 3)
m
出力
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
課題:
関数 installed.packages()
を使用して、現在コンピュータにインストールされているすべてのパッケージの情報を含む
character
行列を作成してください。その内容を確認してください。
R
## create the matrix
ip <- installed.packages()
head(ip)
## try also View(ip)
## number of package
nrow(ip)
## names of all installed packages
rownames(ip)
## type of information we have about each package
colnames(ip)
テストデータとして、大規模なランダムデータ行列を作成することが有用な場合が多くあります。以下の演習では、平均0、標準偏差1の正規分布から抽出したランダムデータを用いて、このような行列を作成することが求められます。これはrnorm()関数で実現可能です。
課題:
正規分布に従うデータ(平均0、標準偏差1)の1000×3の行列を作成してください。
R
set.seed(123)
m <- matrix(rnorm(3000), ncol = 3)
dim(m)
出力
[1] 1000 3
R
head(m)
出力
[,1] [,2] [,3]
[1,] -0.56047565 -0.99579872 -0.5116037
[2,] -0.23017749 -1.03995504 0.2369379
[3,] 1.55870831 -0.01798024 -0.5415892
[4,] 0.07050839 -0.13217513 1.2192276
[5,] 0.12928774 -2.54934277 0.1741359
[6,] 1.71506499 1.04057346 -0.6152683
日付の書式設定
Rの初心者(そして経験豊富な方!)がよく直面する問題の一つに、日付や時刻の情報を、解析で適切かつ利用可能な変数に変換することがあります。
スプレッドシートソフトにおける日付に関する注意点
スプレッドシートにおける日付は、一般的に単一の列に格納されます。これは日付を記録する上で最も自然な方法に思えますが、実際にはベストプラクティスではありません。スプレッドシートアプリケーションは(人間の目には)一見正しく日付を表示しますが、実際の処理や保存方法に問題がある場合があります。日付を「年」「月」「日」の別々の列に、あるいは「年」と「年内日数」の別々の列に格納する方が、はるかに安全なことが多いです。
LibreOffice、Microsoft Excel、OpenOffice、Gnumericなどのスプレッドシートソフトでは、日付のエンコード方式が異なり(多くの場合、互換性がありません)、(同一ソフトであってもバージョンやOSによって異なります)。さらに、Excelでは日付ではないものを日付に変換することがあります( Zeeberg et al. (2004) )。例えば、MAR1、DEC1、OCT4といった名称や識別子などが該当します。そのため、全体的に日付形式を避けることで、こうした問題を特定しやすくなります。
Data Carpentryのレッスンの日付をデータとして扱うセクションでは、スプレッドシートで日付を扱う際の落とし穴について、さらなる知見を提供しています。
lubridate
パッケージ(tidyverse に属します。詳細は
tidyverse のウェブサイト
をご覧ください)の ymd()
関数を使用します。lubridate は
tidyverse
のインストール時に併せてインストールされます。tidyverse
(library(tidyverse))
を読み込むと、コアパッケージ(ほとんどのデータ分析で使用されるパッケージ)が読み込まれます。ただし、lubridate
は tidyverse
のコアパッケージではないため、library(lubridate)
で明示的に読み込む必要があります。
まず、必要なパッケージを読み込みます。
R
library("lubridate")
ymd()は、年・月・日を表すベクトルを受け取り、それをDateベクトルに変換します。Dateは、Rが日付として認識し、それに応じた操作が可能なデータ型です。関数が指定する引数は柔軟ですが、ベストプラクティスとしては、「YYYY-MM-DD」形式の文字列ベクトルを指定します。
日付オブジェクトを作成して、その構造を確認してみましょう。
R
my_date <- ymd("2015-01-01")
str(my_date)
出力
Date[1:1], format: "2015-01-01"
それでは、年、月、日をそれぞれ個別に結合してみましょう。結果は同じになります:
R
# sep indicates the character to use to separate each component
my_date <- ymd(paste("2015", "1", "1", sep = "-"))
str(my_date)
出力
Date[1:1], format: "2015-01-01"
それでは、典型的なデータ操作のパイプラインについて見ていきましょう。以下の小規模なデータセットでは、日付がそれぞれ異なるyear、month、day列に格納されています。
R
x <- data.frame(year = c(1996, 1992, 1987, 1986, 2000, 1990, 2002, 1994, 1997, 1985),
month = c(2, 3, 3, 10, 1, 8, 3, 4, 5, 5),
day = c(24, 8, 1, 5, 8, 17, 13, 10, 11, 24),
value = c(4, 5, 1, 9, 3, 8, 10, 2, 6, 7))
x
出力
year month day value
1 1996 2 24 4
2 1992 3 8 5
3 1987 3 1 1
4 1986 10 5 9
5 2000 1 8 3
6 1990 8 17 8
7 2002 3 13 10
8 1994 4 10 2
9 1997 5 11 6
10 1985 5 24 7
次に、この関数をxデータセットに適用します。まず、paste()を用いてxのyear、month、day列から文字ベクトルを作成します。
R
paste(x$year, x$month, x$day, sep = "-")
出力
[1] "1996-2-24" "1992-3-8" "1987-3-1" "1986-10-5" "2000-1-8" "1990-8-17"
[7] "2002-3-13" "1994-4-10" "1997-5-11" "1985-5-24"
この文字ベクトルは、ymd()の引数として使用できます。
R
ymd(paste(x$year, x$month, x$day, sep = "-"))
出力
[1] "1996-02-24" "1992-03-08" "1987-03-01" "1986-10-05" "2000-01-08"
[6] "1990-08-17" "2002-03-13" "1994-04-10" "1997-05-11" "1985-05-24"
得られたDateベクトルは、xにdateという新しい列として追加できます。
R
x$date <- ymd(paste(x$year, x$month, x$day, sep = "-"))
str(x) # notice the new column, with 'date' as the class
出力
'data.frame': 10 obs. of 5 variables:
$ year : num 1996 1992 1987 1986 2000 ...
$ month: num 2 3 3 10 1 8 3 4 5 5
$ day : num 24 8 1 5 8 17 13 10 11 24
$ value: num 4 5 1 9 3 8 10 2 6 7
$ date : Date, format: "1996-02-24" "1992-03-08" ...
すべて正しく動作したか確認しましょう。新しい列を確認する方法の一つは、summary()を使用することです。
R
summary(x$date)
出力
Min. 1st Qu. Median Mean 3rd Qu. Max.
"1985-05-24" "1988-01-11" "1993-03-24" "1993-03-18" "1997-01-20" "2002-03-13"
なお、ymd()では年、月、日の順で指定する必要があります。例えば、日、月、年の順で指定する場合は、dmy()を使用する必要があります。
R
dmy(paste(x$day, x$month, x$year, sep = "-"))
出力
[1] "1996-02-24" "1992-03-08" "1987-03-01" "1986-10-05" "2000-01-08"
[6] "1990-08-17" "2002-03-13" "1994-04-10" "1997-05-11" "1985-05-24"
lubdridateには、あらゆる日付のバリエーションに対応するための関数が数多く用意されています。
Rオブジェクトの概要
これまでに、次元数や、単一または複数のデータ型を格納できるかどうかが異なる、いくつかの種類のRオブジェクトを見てきました。
-
vector: 1次元(長さがあります)、単一のデータ型です。 -
matrix: 2次元、単一種類のデータ。 -
data.frame: 2次元、1列につき1種類。
リスト
まだ取り上げていませんが、知っておくと便利で、先ほどのまとめから導き出されるデータ型は、リストです。
-
list: 1次元で、各要素は異なるデータ型を持つことができます。
以下に、数値、文字列、行列、データフレーム、および別のリストを含むリストを作成してみましょう。
R
l <- list(1:10, ## numeric
letters, ## character
installed.packages(), ## a matrix
cars, ## a data.frame
list(1, 2, 3)) ## a list
length(l)
出力
[1] 5
R
str(l)
出力
List of 5
$ : int [1:10] 1 2 3 4 5 6 7 8 9 10
$ : chr [1:26] "a" "b" "c" "d" ...
$ : chr [1:261, 1:17] "litedown" "metabolomicsWorkbenchR" "mgcv" "MultiAssayExperiment" ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:261] "litedown" "metabolomicsWorkbenchR" "mgcv" "MultiAssayExperiment" ...
.. ..$ : chr [1:17] "Package" "LibPath" "Version" "Priority" ...
$ :'data.frame': 50 obs. of 2 variables:
..$ speed: num [1:50] 4 4 7 7 8 9 10 10 10 11 ...
..$ dist : num [1:50] 2 10 4 22 16 10 18 26 34 17 ...
$ :List of 3
..$ : num 1
..$ : num 2
..$ : num 3
リストの部分抽出には、新しいサブリストを作成する場合は[]を、そのリストから単一の要素を抽出する場合は[[]](リストに名前が付いている場合はインデックスまたは名前を使用)を使用します。
R
l[[1]] ## first element
出力
[1] 1 2 3 4 5 6 7 8 9 10
R
l[1:2] ## a list of length 2
出力
[[1]]
[1] 1 2 3 4 5 6 7 8 9 10
[[2]]
[1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"
[20] "t" "u" "v" "w" "x" "y" "z"
R
l[1] ## a list of length 1
出力
[[1]]
[1] 1 2 3 4 5 6 7 8 9 10
表形式データのエクスポートと保存
read.table系の関数を用いて、テキスト形式のスプレッドシートをRに読み込む方法を確認しました。data.frameをテキスト形式のスプレッドシートにエクスポートするには、write.table系の関数(write.csv、write.delim、…)を使用します。これらの関数はいずれも、エクスポート対象の変数とエクスポート先のファイルを引数として取ります。例えば、rnaのデータをdata_outputディレクトリ内のmy_rna.csvファイルにエクスポートするには、次のように実行します。
R
write.csv(rna, file = "data_output/my_rna.csv")
この新しいCSVファイルは、Rに慣れていない他の共同研究者とも共有できるようになりました。なお、data.frameの一部のフィールドにはカンマが含まれていますが(例えば「product」列を参照)、Rはデフォルトで各フィールドを引用符で囲むため、列区切りにカンマを使用していても、正しくRに読み込むことができます。
- Rにおける表形式データ
値の長さからもう一方を導き出せるため、行数か列数のどちらか一方があれば十分です。値と行数・列数が整合しない場合にどうなるか試してみてください。↩︎
Content from dplyrを用いたデータの操作と解析
最終更新日:2026-11-07 | ページの編集
所要時間: 150分
概要
質問
- tidyverseメタパッケージを用いたRによるデータ分析
目的
-
dplyrおよびtidyrパッケージの目的を説明してください。 - データ操作に極めて有用な、それらの関数をいくつか説明してください。
- ワイド形式とロング形式のテーブルの概念について解説し、データフレームを一方の形式からもう一方の形式へリシェイプする方法を確認します。
- テーブルの結合方法を説明します。
このエピソードは、Data Carpentriesの「生態学者のためのRによるデータ分析と可視化」レッスンに基づいています。
dplyrとtidyrを用いたデータ操作
角括弧による部分指定は便利ですが、特に複雑な操作では煩雑で読みづらくなることがあります。
データを操作する際、いくつかのパッケージを利用することで作業を大幅に効率化できます。Rのパッケージとは、基本的にはより多くの処理を可能にする追加関数の集合です。これまで使用してきたstr()やdata.frame()などの関数はRに標準で組み込まれています。パッケージをロードすることで、その他の特定の関数を利用できるようになります。パッケージを初めて使用する前に、マシンにインストールしておく必要があります。その後、必要になった際は、それ以降のすべてのRセッションでパッケージをインポートしてください。
パッケージ
dplyrは、データ操作タスクのための強力なツールを提供します。データフレームを直接扱うように設計されており、多くの操作タスクが最適化されています。後述するように、特定の解析や可視化を行うために、データフレームの形状を変更したい場合があります。パッケージ
tidyrは、こうしたデータの形状変更というよくある課題に対応し、データを整然と操作するためのツールを提供します。
ワークショップ終了後にdplyrやtidyrについてさらに詳しく知りたい方は、こちらのdplyrを使った便利なデータ変換チートシートや、こちらのtidyrに関するものをご参照ください。
-
tidyverseパッケージは、tidyr、dplyr、ggplot2、tibbleなど、互いにうまく連携するデータ分析用の便利なパッケージを複数インストールする「アンブレラパッケージ」です。これらのパッケージは、データの操作や解析に役立ちます。これにより、データのサブセット抽出、変換、可視化など、さまざまな処理が可能になります。
セットアップを行ったのであれば、すでにtidyverseパッケージがインストールされているはずです。ライブラリからロードして、インストールされているか確認してください。
R
## load the tidyverse packages, incl. dplyr
library("tidyverse")
エラーメッセージ there is no package called ‘tidyverse’
が表示された場合は、このバージョンの R
向けにパッケージがまだインストールされていません。tidyverse
パッケージをインストールするには、以下を入力してください:
R
BiocManager::install("tidyverse")
tidyverse
パッケージをインストールした場合は、上記の library()
コマンドを使用して、このRセッションでロードするのを忘れないでください!
tidyverseによるデータの読み込み
read.csv()の代わりに、tidyverseパッケージreadrのread_csv()関数(.ではなく_となっている点にご注意ください)を用いてデータを読み込みます。
R
rna <- read_csv("data/rnaseq.csv")
## view the data
rna
出力
# A tibble: 32,428 × 19
gene sample expression organism age sex infection strain time tissue
<chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 Asl GSM254… 1170 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
2 Apod GSM254… 36194 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
3 Cyp2d22 GSM254… 4060 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
4 Klk6 GSM254… 287 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
5 Fcrls GSM254… 85 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
6 Slc2a4 GSM254… 782 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
7 Exd2 GSM254… 1619 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
8 Gjc2 GSM254… 288 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
9 Plp1 GSM254… 43217 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
10 Gnb4 GSM254… 1071 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
# ℹ 32,418 more rows
# ℹ 9 more variables: mouse <dbl>, ENTREZID <dbl>, product <chr>,
# ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
# gene_biotype <chr>, phenotype_description <chr>,
# hsapiens_homolog_associated_gene_name <chr>
データのクラスが現在「tibble」と呼ばれている点にご留意ください。
Tibblesは、前述で導入したデータフレームオブジェクトの挙動の一部を調整したものです。データ構造はデータフレームと非常に似ています。今回の目的においては、相違点は以下の点のみです。
各列の名前の下に、そのデータ型が表示されます。なお、<
dbl>は、小数点を含む数値を保持するように定義されたデータ型です。データの最初の数行と、1画面に収まる分だけの列しか表示されません。
ここからは、最もよく使われるdplyr関数をいくつか学んでいきましょう。
-
select(): 列のサブセット -
filter(): 条件に基づいて行をサブセット化します -
mutate(): 他の列の情報を用いて新しい列を作成します -
group_by()とsummarise(): グループ化されたデータの要約統計量を作成します -
arrange(): 結果を並べ替える -
count(): 離散値をカウントします
列の選択と行のフィルタリング
データフレームの列を選択するには、select()を使用します。この関数の第1引数はデータフレーム(rna)で、それ以降の引数は保持する列です。
R
select(rna, gene, sample, tissue, expression)
出力
# A tibble: 32,428 × 4
gene sample tissue expression
<chr> <chr> <chr> <dbl>
1 Asl GSM2545336 Cerebellum 1170
2 Apod GSM2545336 Cerebellum 36194
3 Cyp2d22 GSM2545336 Cerebellum 4060
4 Klk6 GSM2545336 Cerebellum 287
5 Fcrls GSM2545336 Cerebellum 85
6 Slc2a4 GSM2545336 Cerebellum 782
7 Exd2 GSM2545336 Cerebellum 1619
8 Gjc2 GSM2545336 Cerebellum 288
9 Plp1 GSM2545336 Cerebellum 43217
10 Gnb4 GSM2545336 Cerebellum 1071
# ℹ 32,418 more rows
特定の列以外をすべて選択するには、除外したい変数の前に「-」を付けます。
R
select(rna, -tissue, -organism)
出力
# A tibble: 32,428 × 17
gene sample expression age sex infection strain time mouse ENTREZID
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl>
1 Asl GSM2545… 1170 8 Fema… Influenz… C57BL… 8 14 109900
2 Apod GSM2545… 36194 8 Fema… Influenz… C57BL… 8 14 11815
3 Cyp2d22 GSM2545… 4060 8 Fema… Influenz… C57BL… 8 14 56448
4 Klk6 GSM2545… 287 8 Fema… Influenz… C57BL… 8 14 19144
5 Fcrls GSM2545… 85 8 Fema… Influenz… C57BL… 8 14 80891
6 Slc2a4 GSM2545… 782 8 Fema… Influenz… C57BL… 8 14 20528
7 Exd2 GSM2545… 1619 8 Fema… Influenz… C57BL… 8 14 97827
8 Gjc2 GSM2545… 288 8 Fema… Influenz… C57BL… 8 14 118454
9 Plp1 GSM2545… 43217 8 Fema… Influenz… C57BL… 8 14 18823
10 Gnb4 GSM2545… 1071 8 Fema… Influenz… C57BL… 8 14 14696
# ℹ 32,418 more rows
# ℹ 7 more variables: product <chr>, ensembl_gene_id <chr>,
# external_synonym <chr>, chromosome_name <chr>, gene_biotype <chr>,
# phenotype_description <chr>, hsapiens_homolog_associated_gene_name <chr>
これにより、rna 内の tissue と
organism を除くすべての変数が選択されます。
特定の条件に基づいて行を選択するには、filter()を使用します。
R
filter(rna, sex == "Male")
出力
# A tibble: 14,740 × 19
gene sample expression organism age sex infection strain time tissue
<chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 Asl GSM254… 626 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
2 Apod GSM254… 13021 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
3 Cyp2d22 GSM254… 2171 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
4 Klk6 GSM254… 448 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
5 Fcrls GSM254… 180 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
6 Slc2a4 GSM254… 313 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
7 Exd2 GSM254… 2366 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
8 Gjc2 GSM254… 310 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
9 Plp1 GSM254… 53126 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
10 Gnb4 GSM254… 1355 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
# ℹ 14,730 more rows
# ℹ 9 more variables: mouse <dbl>, ENTREZID <dbl>, product <chr>,
# ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
# gene_biotype <chr>, phenotype_description <chr>,
# hsapiens_homolog_associated_gene_name <chr>
R
filter(rna, sex == "Male" & infection == "NonInfected")
出力
# A tibble: 4,422 × 19
gene sample expression organism age sex infection strain time tissue
<chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 Asl GSM254… 535 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
2 Apod GSM254… 13668 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
3 Cyp2d22 GSM254… 2008 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
4 Klk6 GSM254… 1101 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
5 Fcrls GSM254… 375 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
6 Slc2a4 GSM254… 249 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
7 Exd2 GSM254… 3126 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
8 Gjc2 GSM254… 791 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
9 Plp1 GSM254… 98658 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
10 Gnb4 GSM254… 2437 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
# ℹ 4,412 more rows
# ℹ 9 more variables: mouse <dbl>, ENTREZID <dbl>, product <chr>,
# ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
# gene_biotype <chr>, phenotype_description <chr>,
# hsapiens_homolog_associated_gene_name <chr>
それでは、このデータセットで解析されたマウス遺伝子のヒト相同遺伝子に関心があるとしましょう。この情報は、rna
ティブルの最後の列である
hsapiens_homolog_associated_gene_name
に格納されています。これを分かりやすく可視化するため、gene
と hsapiens_homolog_associated_gene_name
の2つの列のみを含む新しいテーブルを作成します。
R
genes <- select(rna, gene, hsapiens_homolog_associated_gene_name)
genes
出力
# A tibble: 32,428 × 2
gene hsapiens_homolog_associated_gene_name
<chr> <chr>
1 Asl ASL
2 Apod APOD
3 Cyp2d22 CYP2D6
4 Klk6 KLK6
5 Fcrls FCRL2
6 Slc2a4 SLC2A4
7 Exd2 EXD2
8 Gjc2 GJC2
9 Plp1 PLP1
10 Gnb4 GNB4
# ℹ 32,418 more rows
マウスの遺伝子の中には、ヒトの相同遺伝子が存在しないものがあります。これらは、filter()
および、あるものが NA であるかどうかを判定する
is.na() 関数を用いて取得できます。
R
filter(genes, is.na(hsapiens_homolog_associated_gene_name))
出力
# A tibble: 4,290 × 2
gene hsapiens_homolog_associated_gene_name
<chr> <chr>
1 Prodh <NA>
2 Tssk5 <NA>
3 Vmn2r1 <NA>
4 Gm10654 <NA>
5 Hexa <NA>
6 Sult1a1 <NA>
7 Gm6277 <NA>
8 Tmem198b <NA>
9 Adam1a <NA>
10 Ebp <NA>
# ℹ 4,280 more rows
ヒトの相同遺伝子を持つマウス遺伝子のみを保持したい場合は、結果を否定する「!」記号を挿入することで、hsapiens_homolog_associated_gene_nameがNAではないすべての行を抽出できます。
R
filter(genes, !is.na(hsapiens_homolog_associated_gene_name))
出力
# A tibble: 28,138 × 2
gene hsapiens_homolog_associated_gene_name
<chr> <chr>
1 Asl ASL
2 Apod APOD
3 Cyp2d22 CYP2D6
4 Klk6 KLK6
5 Fcrls FCRL2
6 Slc2a4 SLC2A4
7 Exd2 EXD2
8 Gjc2 GJC2
9 Plp1 PLP1
10 Gnb4 GNB4
# ℹ 28,128 more rows
パイプ
選択とフィルタリングを同時に行いたい場合はどうすればよいでしょうか。その方法には、中間ステップを使う、ネストした関数を使う、あるいはパイプを使うという3通りがあります。
中間ステップを挟むことで、一時的なデータフレームを作成し、それを次の関数への入力として次のように利用します。
R
rna2 <- filter(rna, sex == "Male")
rna3 <- select(rna2, gene, sample, tissue, expression)
rna3
出力
# A tibble: 14,740 × 4
gene sample tissue expression
<chr> <chr> <chr> <dbl>
1 Asl GSM2545340 Cerebellum 626
2 Apod GSM2545340 Cerebellum 13021
3 Cyp2d22 GSM2545340 Cerebellum 2171
4 Klk6 GSM2545340 Cerebellum 448
5 Fcrls GSM2545340 Cerebellum 180
6 Slc2a4 GSM2545340 Cerebellum 313
7 Exd2 GSM2545340 Cerebellum 2366
8 Gjc2 GSM2545340 Cerebellum 310
9 Plp1 GSM2545340 Cerebellum 53126
10 Gnb4 GSM2545340 Cerebellum 1355
# ℹ 14,730 more rows
これは読みやすいものの、個別に名前を付ける必要がある中間オブジェクトが大量に発生し、作業スペースが煩雑になることがあります。手順が複数に及ぶと、把握しにくくなる場合があります。
また、次のように関数をネスト(つまり、ある関数の中に別の関数を入れる)することもできます。
R
rna3 <- select(filter(rna, sex == "Male"), gene, sample, tissue, expression)
rna3
出力
# A tibble: 14,740 × 4
gene sample tissue expression
<chr> <chr> <chr> <dbl>
1 Asl GSM2545340 Cerebellum 626
2 Apod GSM2545340 Cerebellum 13021
3 Cyp2d22 GSM2545340 Cerebellum 2171
4 Klk6 GSM2545340 Cerebellum 448
5 Fcrls GSM2545340 Cerebellum 180
6 Slc2a4 GSM2545340 Cerebellum 313
7 Exd2 GSM2545340 Cerebellum 2366
8 Gjc2 GSM2545340 Cerebellum 310
9 Plp1 GSM2545340 Cerebellum 53126
10 Gnb4 GSM2545340 Cerebellum 1355
# ℹ 14,730 more rows
これは便利ですが、関数が入れ子になりすぎると読みにくくなることがあります。Rは式を内側から外側に向かって評価するためです(この場合は、フィルタリングしてから選択を行います)。
最後のオプションであるパイプは、Rに最近追加された機能です。パイプを使うと、ある関数の出力をそのまま次の関数に送ることができるため、同じデータセットに対して多くの処理を行う際に便利です。
Rにおけるパイプは
%>%(magrittr
パッケージで利用可能)や |>(base
R経由)のような形式になります。RStudioを使用している場合、PCであれば
Ctrl + Shift + M、Macであれば
Cmd + Shift + M と入力できます。
上記のコードでは、パイプを用いて、まずrnaデータセットをfilter()に通してsexがMaleの行を残し、次にselect()に通してgene、sample、tissue、expressionの列のみを残しています。
パイプ |>
は、左側のオブジェクトを受け取り、それを右側の関数の第1引数としてそのまま渡します。そのため、filter()
や select()
関数の引数としてデータフレームを明示的に指定する必要はもうありません。
R
rna |>
filter(sex == "Male") |>
select(gene, sample, tissue, expression)
出力
# A tibble: 14,740 × 4
gene sample tissue expression
<chr> <chr> <chr> <dbl>
1 Asl GSM2545340 Cerebellum 626
2 Apod GSM2545340 Cerebellum 13021
3 Cyp2d22 GSM2545340 Cerebellum 2171
4 Klk6 GSM2545340 Cerebellum 448
5 Fcrls GSM2545340 Cerebellum 180
6 Slc2a4 GSM2545340 Cerebellum 313
7 Exd2 GSM2545340 Cerebellum 2366
8 Gjc2 GSM2545340 Cerebellum 310
9 Plp1 GSM2545340 Cerebellum 53126
10 Gnb4 GSM2545340 Cerebellum 1355
# ℹ 14,730 more rows
パイプを「then」という単語のように読むと分かりやすいと感じる方もいるかもしれません。例えば、上記の例では、データフレーム
rna を取り出し、then sex == "Male"
を持つ行を filter し、then 列
gene、sample、tissue、expression
を select しました。
dplyr
関数単体ではいたってシンプルですが、パイプを用いて線形なワークフローに組み合わせることで、データフレームに対してより複雑な操作を行うことができます。
この小規模版のデータを用いて新しいオブジェクトを作成したい場合は、それに新しい名前を付けることができます。
R
rna3 <- rna |>
filter(sex == "Male") |>
select(gene, sample, tissue, expression)
rna3
出力
# A tibble: 14,740 × 4
gene sample tissue expression
<chr> <chr> <chr> <dbl>
1 Asl GSM2545340 Cerebellum 626
2 Apod GSM2545340 Cerebellum 13021
3 Cyp2d22 GSM2545340 Cerebellum 2171
4 Klk6 GSM2545340 Cerebellum 448
5 Fcrls GSM2545340 Cerebellum 180
6 Slc2a4 GSM2545340 Cerebellum 313
7 Exd2 GSM2545340 Cerebellum 2366
8 Gjc2 GSM2545340 Cerebellum 310
9 Plp1 GSM2545340 Cerebellum 53126
10 Gnb4 GSM2545340 Cerebellum 1355
# ℹ 14,730 more rows
課題:
パイプを用いてrnaデータをサブセットし、時刻0における雌マウスの観測値のうち、遺伝子発現量が50000を超えるものを抽出し、gene、sample、time、expression、ageの各列のみを保持してください。
R
rna |>
filter(expression > 50000,
sex == "Female",
time == 0 ) |>
select(gene, sample, time, expression, age)
出力
# A tibble: 9 × 5
gene sample time expression age
<chr> <chr> <dbl> <dbl> <dbl>
1 Plp1 GSM2545337 0 101241 8
2 Atp1b1 GSM2545337 0 53260 8
3 Plp1 GSM2545338 0 96534 8
4 Atp1b1 GSM2545338 0 50614 8
5 Plp1 GSM2545348 0 102790 8
6 Atp1b1 GSM2545348 0 59544 8
7 Plp1 GSM2545353 0 71237 8
8 Glul GSM2545353 0 52451 8
9 Atp1b1 GSM2545353 0 61451 8
列の作成・変更(mutate)
既存の列の値に基づいて、例えば単位換算を行ったり、2つの列の値の比率を求めたりするために、新しい列を作成したいことがよくあります。ここではmutate()を使用します。
時間(単位:時間)の新しい列を作成するには:
R
rna |>
mutate(time_hours = time * 24) |>
select(time, time_hours)
出力
# A tibble: 32,428 × 2
time time_hours
<dbl> <dbl>
1 8 192
2 8 192
3 8 192
4 8 192
5 8 192
6 8 192
7 8 192
8 8 192
9 8 192
10 8 192
# ℹ 32,418 more rows
また、mutate()
の同一呼び出し内で、最初の新規列に基づいて2つ目の新規列を作成することもできます。
R
rna |>
mutate(time_hours = time * 24,
time_mn = time_hours * 60) |>
select(time, time_hours, time_mn)
出力
# A tibble: 32,428 × 3
time time_hours time_mn
<dbl> <dbl> <dbl>
1 8 192 11520
2 8 192 11520
3 8 192 11520
4 8 192 11520
5 8 192 11520
6 8 192 11520
7 8 192 11520
8 8 192 11520
9 8 192 11520
10 8 192 11520
# ℹ 32,418 more rows
課題
rnaのデータから、以下の条件を満たす新しいデータフレームを作成してください。gene、chromosome_name、phenotype_description、sample、およびexpressionの列のみを含むものとします。発現値は対数変換してください。このデータフレームには、性染色体上に位置し、phenotype_descriptionに関連付けられており、かつ対数発現値が5を超える遺伝子のみを含める必要があります。
ヒント: このデータフレームを生成するには、コマンドをどのような順序で並べればよいか考えてみてください!
R
rna |>
mutate(expression = log(expression)) |>
select(gene, chromosome_name, phenotype_description, sample, expression) |>
filter(chromosome_name == "X" | chromosome_name == "Y") |>
filter(!is.na(phenotype_description)) |>
filter(expression > 5)
出力
# A tibble: 649 × 5
gene chromosome_name phenotype_description sample expression
<chr> <chr> <chr> <chr> <dbl>
1 Plp1 X abnormal CNS glial cell morphology GSM25… 10.7
2 Slc7a3 X decreased body length GSM25… 5.46
3 Plxnb3 X abnormal coat appearance GSM25… 6.58
4 Rbm3 X abnormal liver morphology GSM25… 9.32
5 Cfp X abnormal cardiovascular system phys… GSM25… 6.18
6 Ebp X abnormal embryonic erythrocyte morp… GSM25… 6.68
7 Cd99l2 X abnormal cellular extravasation GSM25… 8.04
8 Piga X abnormal brain development GSM25… 6.06
9 Pim2 X decreased T cell proliferation GSM25… 7.11
10 Itm2a X no abnormal phenotype detected GSM25… 7.48
# ℹ 639 more rows
分割・適用・結合によるデータ分析
多くのデータ分析タスクは、split-apply-combine
パラダイムを用いて取り組むことができます。つまり、データをグループに分割し、各グループに何らかの分析を適用し、その結果を統合するという手法です。dplyr
では、group_by()
関数を用いることで、これを非常に簡単に実現できます。
R
rna |>
group_by(gene)
出力
# A tibble: 32,428 × 19
# Groups: gene [1,474]
gene sample expression organism age sex infection strain time tissue
<chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 Asl GSM254… 1170 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
2 Apod GSM254… 36194 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
3 Cyp2d22 GSM254… 4060 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
4 Klk6 GSM254… 287 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
5 Fcrls GSM254… 85 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
6 Slc2a4 GSM254… 782 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
7 Exd2 GSM254… 1619 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
8 Gjc2 GSM254… 288 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
9 Plp1 GSM254… 43217 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
10 Gnb4 GSM254… 1071 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
# ℹ 32,418 more rows
# ℹ 9 more variables: mouse <dbl>, ENTREZID <dbl>, product <chr>,
# ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
# gene_biotype <chr>, phenotype_description <chr>,
# hsapiens_homolog_associated_gene_name <chr>
group_by()関数はデータ処理を行わず、データをサブセットにグループ化します。上記の例では、初期のtibbleに含まれる32428個の観測値が、gene変数に基づいて1474個のグループに分割されます。
同様に、tibbleをサンプルごとにグループ化することもできます。
R
rna |>
group_by(sample)
出力
# A tibble: 32,428 × 19
# Groups: sample [22]
gene sample expression organism age sex infection strain time tissue
<chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 Asl GSM254… 1170 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
2 Apod GSM254… 36194 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
3 Cyp2d22 GSM254… 4060 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
4 Klk6 GSM254… 287 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
5 Fcrls GSM254… 85 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
6 Slc2a4 GSM254… 782 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
7 Exd2 GSM254… 1619 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
8 Gjc2 GSM254… 288 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
9 Plp1 GSM254… 43217 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
10 Gnb4 GSM254… 1071 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
# ℹ 32,418 more rows
# ℹ 9 more variables: mouse <dbl>, ENTREZID <dbl>, product <chr>,
# ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
# gene_biotype <chr>, phenotype_description <chr>,
# hsapiens_homolog_associated_gene_name <chr>
ここでは、32428 個の観測値からなる初期の tibble
を、sample 変数に基づいて 22 個のグループに分割します。
データをグループ化した後は、以降の処理は各グループに対して個別に適用されます。
summarise()関数
group_by()はsummarise()と併用されることが多く、これにより各グループがそのグループの1行の要約にまとめられます。
group_by()は、要約統計量を算出したいカテゴリカル変数を含む列名を引数として受け取ります。したがって、遺伝子ごとの平均値
expression を算出するには:
R
rna |>
group_by(gene) |>
summarise(mean_expression = mean(expression))
出力
# A tibble: 1,474 × 2
gene mean_expression
<chr> <dbl>
1 AI504432 1053.
2 AW046200 131.
3 AW551984 295.
4 Aamp 4751.
5 Abca12 4.55
6 Abcc8 2498.
7 Abhd14a 525.
8 Abi2 4909.
9 Abi3bp 1002.
10 Abl2 2124.
# ℹ 1,464 more rows
また、各サンプルにおける全遺伝子の平均発現量を算出したい場合もあるでしょう。
R
rna |>
group_by(sample) |>
summarise(mean_expression = mean(expression))
出力
# A tibble: 22 × 2
sample mean_expression
<chr> <dbl>
1 GSM2545336 2062.
2 GSM2545337 1766.
3 GSM2545338 1668.
4 GSM2545339 1696.
5 GSM2545340 1682.
6 GSM2545341 1638.
7 GSM2545342 1594.
8 GSM2545343 2107.
9 GSM2545344 1712.
10 GSM2545345 1700.
# ℹ 12 more rows
また、複数の列でグループ化することも可能です。
R
rna |>
group_by(gene, infection, time) |>
summarise(mean_expression = mean(expression))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene, infection, and time.
ℹ Output is grouped by gene and infection.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, infection, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 4,422 × 4
# Groups: gene, infection [2,948]
gene infection time mean_expression
<chr> <chr> <dbl> <dbl>
1 AI504432 InfluenzaA 4 1104.
2 AI504432 InfluenzaA 8 1014
3 AI504432 NonInfected 0 1034.
4 AW046200 InfluenzaA 4 152.
5 AW046200 InfluenzaA 8 81
6 AW046200 NonInfected 0 155.
7 AW551984 InfluenzaA 4 302.
8 AW551984 InfluenzaA 8 342.
9 AW551984 NonInfected 0 238
10 Aamp InfluenzaA 4 4870
# ℹ 4,412 more rows
データをグループ化した後は、複数の変数を同時に(必ずしも同一の変数である必要はありません)集計することも可能です。例えば、遺伝子別および条件別に中央値
expression を示す列を追加できます。
R
rna |>
group_by(gene, infection, time) |>
summarise(mean_expression = mean(expression),
median_expression = median(expression))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene, infection, and time.
ℹ Output is grouped by gene and infection.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, infection, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 4,422 × 5
# Groups: gene, infection [2,948]
gene infection time mean_expression median_expression
<chr> <chr> <dbl> <dbl> <dbl>
1 AI504432 InfluenzaA 4 1104. 1094.
2 AI504432 InfluenzaA 8 1014 985
3 AI504432 NonInfected 0 1034. 1016
4 AW046200 InfluenzaA 4 152. 144.
5 AW046200 InfluenzaA 8 81 82
6 AW046200 NonInfected 0 155. 163
7 AW551984 InfluenzaA 4 302. 245
8 AW551984 InfluenzaA 8 342. 287
9 AW551984 NonInfected 0 238 265
10 Aamp InfluenzaA 4 4870 4708
# ℹ 4,412 more rows
課題
時点ごとの遺伝子「Dok3」の平均発現量を算出してください。
R
rna |>
filter(gene == "Dok3") |>
group_by(time) |>
summarise(mean = mean(expression))
出力
# A tibble: 3 × 2
time mean
<dbl> <dbl>
1 0 169
2 4 156.
3 8 61
カウント
データを扱う際、各因子または因子の組み合わせごとに観測数を確認したいことがよくあります。この操作には、dplyr
で count()
を提供します。例えば、感染サンプルと非感染サンプルそれぞれのデータ行数を数えたい場合は、次のように記述します。
R
rna |>
count(infection)
出力
# A tibble: 2 × 2
infection n
<chr> <int>
1 InfluenzaA 22110
2 NonInfected 10318
count()
関数は、すでに見てきた手法の簡略表記です。具体的には、変数でグループ化し、そのグループ内の観測数を集計するものです。言い換えれば、rna |> count(infection)
は以下と同等です。
R
rna |>
group_by(infection) |>
summarise(n = n())
出力
# A tibble: 2 × 2
infection n
<chr> <int>
1 InfluenzaA 22110
2 NonInfected 10318
前述の例では、count()
を用いて1つの因子(すなわち
infection)の行数/観測数をカウントしています。infection
と time
のような因子の組み合わせをカウントしたい場合は、count()
の引数として第1因子と第2因子を指定します:
R
rna |>
count(infection, time)
出力
# A tibble: 3 × 3
infection time n
<chr> <dbl> <int>
1 InfluenzaA 4 11792
2 InfluenzaA 8 10318
3 NonInfected 0 10318
これは以下と同等です。
R
rna |>
group_by(infection, time) |>
summarise(n = n())
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by infection and time.
ℹ Output is grouped by infection.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(infection, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 3 × 3
# Groups: infection [2]
infection time n
<chr> <dbl> <int>
1 InfluenzaA 4 11792
2 InfluenzaA 8 10318
3 NonInfected 0 10318
比較を容易にするために、結果を並べ替えると便利な場合があります。テーブルの並べ替えには
arrange()
を使用できます。例えば、上記のテーブルを時間順に並べ替えたい場合などが考えられます。
R
rna |>
count(infection, time) |>
arrange(time)
出力
# A tibble: 3 × 3
infection time n
<chr> <dbl> <int>
1 NonInfected 0 10318
2 InfluenzaA 4 11792
3 InfluenzaA 8 10318
または件数で:
R
rna |>
count(infection, time) |>
arrange(n)
出力
# A tibble: 3 × 3
infection time n
<chr> <dbl> <int>
1 InfluenzaA 8 10318
2 NonInfected 0 10318
3 InfluenzaA 4 11792
降順に並べ替えるには、desc()関数を追加する必要があります。
R
rna |>
count(infection, time) |>
arrange(desc(n))
出力
# A tibble: 3 × 3
infection time n
<chr> <dbl> <int>
1 InfluenzaA 4 11792
2 InfluenzaA 8 10318
3 NonInfected 0 10318
課題
- 各サンプルで解析された遺伝子数はいくつですか?
-
group_by()とsummarise()を用いて、各サンプルのシーケンシング深度(全カウントの合計)を算出してください。最もシーケンシング深度が高いサンプルはどれですか? - サンプルを1つ選び、バイオタイプごとの遺伝子数を算出してください。
- 「abnormal DNA methylation」という表現型の説明に関連する遺伝子を特定し、時刻0、4、8におけるそれらの平均発現量(対数スケール)を算出してください。
R
## 1.
rna |>
count(sample)
出力
# A tibble: 22 × 2
sample n
<chr> <int>
1 GSM2545336 1474
2 GSM2545337 1474
3 GSM2545338 1474
4 GSM2545339 1474
5 GSM2545340 1474
6 GSM2545341 1474
7 GSM2545342 1474
8 GSM2545343 1474
9 GSM2545344 1474
10 GSM2545345 1474
# ℹ 12 more rows
R
## 2.
rna |>
group_by(sample) |>
summarise(seq_depth = sum(expression)) |>
arrange(desc(seq_depth))
出力
# A tibble: 22 × 2
sample seq_depth
<chr> <dbl>
1 GSM2545350 3255566
2 GSM2545352 3216163
3 GSM2545343 3105652
4 GSM2545336 3039671
5 GSM2545380 3036098
6 GSM2545353 2953249
7 GSM2545348 2913678
8 GSM2545362 2913517
9 GSM2545351 2782464
10 GSM2545349 2758006
# ℹ 12 more rows
R
## 3.
rna |>
filter(sample == "GSM2545336") |>
count(gene_biotype) |>
arrange(desc(n))
出力
# A tibble: 13 × 2
gene_biotype n
<chr> <int>
1 protein_coding 1321
2 lncRNA 69
3 processed_pseudogene 59
4 miRNA 7
5 snoRNA 5
6 TEC 4
7 polymorphic_pseudogene 2
8 unprocessed_pseudogene 2
9 IG_C_gene 1
10 scaRNA 1
11 transcribed_processed_pseudogene 1
12 transcribed_unitary_pseudogene 1
13 transcribed_unprocessed_pseudogene 1
R
## 4.
rna |>
filter(phenotype_description == "abnormal DNA methylation") |>
group_by(gene, time) |>
summarise(mean_expression = mean(log(expression))) |>
arrange()
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 6 × 3
# Groups: gene [2]
gene time mean_expression
<chr> <dbl> <dbl>
1 Xist 0 6.95
2 Xist 4 6.34
3 Xist 8 7.13
4 Zdbf2 0 6.27
5 Zdbf2 4 6.27
6 Zdbf2 8 6.19
データの再整形
rna ティブルの行には、他の2つの変数 gene と
sample
の組み合わせに関連付けられた発現値(単位)が含まれています。
その他の列はすべて、サンプル(生物種、年齢、性別、…)または遺伝子(gene_biotype、ENTREZ_ID、産物、…)を表す変数に対応しています。遺伝子やサンプルによって変化しない変数は、すべての行で同一の値となります。
R
rna |>
arrange(gene)
出力
# A tibble: 32,428 × 19
gene sample expression organism age sex infection strain time tissue
<chr> <chr> <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 AI504432 GSM25… 1230 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
2 AI504432 GSM25… 1085 Mus mus… 8 Fema… NonInfec… C57BL… 0 Cereb…
3 AI504432 GSM25… 969 Mus mus… 8 Fema… NonInfec… C57BL… 0 Cereb…
4 AI504432 GSM25… 1284 Mus mus… 8 Fema… Influenz… C57BL… 4 Cereb…
5 AI504432 GSM25… 966 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
6 AI504432 GSM25… 918 Mus mus… 8 Male Influenz… C57BL… 8 Cereb…
7 AI504432 GSM25… 985 Mus mus… 8 Fema… Influenz… C57BL… 8 Cereb…
8 AI504432 GSM25… 972 Mus mus… 8 Male NonInfec… C57BL… 0 Cereb…
9 AI504432 GSM25… 1000 Mus mus… 8 Fema… Influenz… C57BL… 4 Cereb…
10 AI504432 GSM25… 816 Mus mus… 8 Male Influenz… C57BL… 4 Cereb…
# ℹ 32,418 more rows
# ℹ 9 more variables: mouse <dbl>, ENTREZID <dbl>, product <chr>,
# ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
# gene_biotype <chr>, phenotype_description <chr>,
# hsapiens_homolog_associated_gene_name <chr>
この構造は、1つの列にすべての値が含まれ、他の列にその値のコンテキストが列挙されているため、long-formatと呼ばれます。
場合によっては、long-formatは必ずしも「人間が読みやすい」とは言えず、データをよりコンパクトに表現できる別の形式であるwide-formatが好まれます。これは通常、科学者が行列形式(行が遺伝子、列がサンプルを表す形式)で見慣れている遺伝子発現値の場合に当てはまります。
この形式であれば、サンプル内およびサンプル間における遺伝子発現レベルの関係を容易に解析できるようになります。
出力
# A tibble: 1,474 × 23
gene GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Asl 1170 361 400 586 626 988
2 Apod 36194 10347 9173 10620 13021 29594
3 Cyp2d22 4060 1616 1603 1901 2171 3349
4 Klk6 287 629 641 578 448 195
5 Fcrls 85 233 244 237 180 38
6 Slc2a4 782 231 248 265 313 786
7 Exd2 1619 2288 2235 2513 2366 1359
8 Gjc2 288 595 568 551 310 146
9 Plp1 43217 101241 96534 58354 53126 27173
10 Gnb4 1071 1791 1867 1430 1355 798
# ℹ 1,464 more rows
# ℹ 16 more variables: GSM2545342 <dbl>, GSM2545343 <dbl>, GSM2545344 <dbl>,
# GSM2545345 <dbl>, GSM2545346 <dbl>, GSM2545347 <dbl>, GSM2545348 <dbl>,
# GSM2545349 <dbl>, GSM2545350 <dbl>, GSM2545351 <dbl>, GSM2545352 <dbl>,
# GSM2545353 <dbl>, GSM2545354 <dbl>, GSM2545362 <dbl>, GSM2545363 <dbl>,
# GSM2545380 <dbl>
遺伝子発現値をrnaからワイド形式に変換するには、sample列の値を列変数名とする新しいテーブルを作成する必要があります。
ここでのポイントは、引き続き整然データ構造を採用しつつ、関心のある観測項目に合わせてデータを再構成した点です。具体的には、遺伝子ごと・サンプルごとに記録するのではなく、遺伝子ごとの発現レベルとしてまとめています。
逆の変換は、列名を新しい変数の値に変換することです。
これら2つの変換は、tidyr 関数の
pivot_longer() と pivot_wider()
の2つで行えます(詳細は tidyr
ウェブサイトのピボットに関する記事 を参照してください)。
データをより幅広の形式にピボットする
rnaの先頭3列を選択し、pivot_wider()を使ってデータをワイド形式に変換しましょう。
R
rna_exp <- rna |>
select(gene, sample, expression)
rna_exp
出力
# A tibble: 32,428 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Apod GSM2545336 36194
3 Cyp2d22 GSM2545336 4060
4 Klk6 GSM2545336 287
5 Fcrls GSM2545336 85
6 Slc2a4 GSM2545336 782
7 Exd2 GSM2545336 1619
8 Gjc2 GSM2545336 288
9 Plp1 GSM2545336 43217
10 Gnb4 GSM2545336 1071
# ℹ 32,418 more rows
pivot_widerは、3つの主要な引数を取ります。
- 変換対象のデータ;
-
names_from: 値が新しい列名となる列; -
values_from: 新しい列に値が代入される列です。

rna data.
R
rna_wide <- rna_exp |>
pivot_wider(names_from = sample,
values_from = expression)
rna_wide
出力
# A tibble: 1,474 × 23
gene GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Asl 1170 361 400 586 626 988
2 Apod 36194 10347 9173 10620 13021 29594
3 Cyp2d22 4060 1616 1603 1901 2171 3349
4 Klk6 287 629 641 578 448 195
5 Fcrls 85 233 244 237 180 38
6 Slc2a4 782 231 248 265 313 786
7 Exd2 1619 2288 2235 2513 2366 1359
8 Gjc2 288 595 568 551 310 146
9 Plp1 43217 101241 96534 58354 53126 27173
10 Gnb4 1071 1791 1867 1430 1355 798
# ℹ 1,464 more rows
# ℹ 16 more variables: GSM2545342 <dbl>, GSM2545343 <dbl>, GSM2545344 <dbl>,
# GSM2545345 <dbl>, GSM2545346 <dbl>, GSM2545347 <dbl>, GSM2545348 <dbl>,
# GSM2545349 <dbl>, GSM2545350 <dbl>, GSM2545351 <dbl>, GSM2545352 <dbl>,
# GSM2545353 <dbl>, GSM2545354 <dbl>, GSM2545362 <dbl>, GSM2545363 <dbl>,
# GSM2545380 <dbl>
なお、デフォルトでは、pivot_wider() 関数は欠損値に対して
NA を追加します。
何らかの理由で、特定のサンプルにおいて一部の遺伝子の発現値が欠損している状況を想定してみましょう。以下の架空の例では、遺伝子Cyp2d22の発現値は、サンプルGSM2545338において1つしかありません。
R
rna_with_missing_values <- rna |>
select(gene, sample, expression) |>
filter(gene %in% c("Asl", "Apod", "Cyp2d22")) |>
filter(sample %in% c("GSM2545336", "GSM2545337", "GSM2545338")) |>
arrange(sample) |>
filter(!(gene == "Cyp2d22" & sample != "GSM2545338"))
rna_with_missing_values
出力
# A tibble: 7 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Apod GSM2545336 36194
3 Asl GSM2545337 361
4 Apod GSM2545337 10347
5 Asl GSM2545338 400
6 Apod GSM2545338 9173
7 Cyp2d22 GSM2545338 1603
デフォルトでは、pivot_wider()関数は欠損値に対してNAを追加します。これは、pivot_wider()関数のvalues_fill引数で指定可能です。
R
rna_with_missing_values |>
pivot_wider(names_from = sample,
values_from = expression)
出力
# A tibble: 3 × 4
gene GSM2545336 GSM2545337 GSM2545338
<chr> <dbl> <dbl> <dbl>
1 Asl 1170 361 400
2 Apod 36194 10347 9173
3 Cyp2d22 NA NA 1603
R
rna_with_missing_values |>
pivot_wider(names_from = sample,
values_from = expression,
values_fill = 0)
出力
# A tibble: 3 × 4
gene GSM2545336 GSM2545337 GSM2545338
<chr> <dbl> <dbl> <dbl>
1 Asl 1170 361 400
2 Apod 36194 10347 9173
3 Cyp2d22 0 0 1603
データを縦長形式にピボットする
逆に、列名を利用して、それらを一対の新しい変数に変換する場合もあります。一方の変数は列名を値として表し、もう一方の変数にはそれまで列名に関連付けられていた値が格納されます。
pivot_longer()は、4つの主要な引数を取ります。
- 変換対象のデータ;
-
names_to: 新たに作成し、現在の列名を格納したい新しい列名です。 -
values_to: 新たに作成し、現在の値を格納したい列名です。 -
names_toおよびvalues_to変数に格納(または削除)する列名です。

rna data.
rna_wideからrna_longを再作成するには、sampleというキーとexpressionという値を作成し、gene以外のすべての列をキー変数として使用します。ここでは、gene列をマイナス記号付きで除外します。
ここでは、新しい変数名をクォートで囲む点にご留意ください。
R
rna_long <- rna_wide |>
pivot_longer(names_to = "sample",
values_to = "expression",
-gene)
rna_long
出力
# A tibble: 32,428 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Asl GSM2545337 361
3 Asl GSM2545338 400
4 Asl GSM2545339 586
5 Asl GSM2545340 626
6 Asl GSM2545341 988
7 Asl GSM2545342 836
8 Asl GSM2545343 535
9 Asl GSM2545344 586
10 Asl GSM2545345 597
# ℹ 32,418 more rows
また、含める列を指定する仕様を用いることもできました。識別用の列が多数ある場合に有用で、除外する列よりも取得する列を指定する方が簡単です。ここでは、starts_with()関数を使えば、すべてを列挙することなくサンプル名を取得できます!別の方法としては、:演算子を使うことも考えられます!
R
rna_wide |>
pivot_longer(names_to = "sample",
values_to = "expression",
cols = starts_with("GSM"))
出力
# A tibble: 32,428 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Asl GSM2545337 361
3 Asl GSM2545338 400
4 Asl GSM2545339 586
5 Asl GSM2545340 626
6 Asl GSM2545341 988
7 Asl GSM2545342 836
8 Asl GSM2545343 535
9 Asl GSM2545344 586
10 Asl GSM2545345 597
# ℹ 32,418 more rows
R
rna_wide |>
pivot_longer(names_to = "sample",
values_to = "expression",
GSM2545336:GSM2545380)
出力
# A tibble: 32,428 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Asl GSM2545337 361
3 Asl GSM2545338 400
4 Asl GSM2545339 586
5 Asl GSM2545340 626
6 Asl GSM2545341 988
7 Asl GSM2545342 836
8 Asl GSM2545343 535
9 Asl GSM2545344 586
10 Asl GSM2545345 597
# ℹ 32,418 more rows
なお、ワイド形式で欠損値があった場合、NA
は新しいロング形式にも含まれる点にご注意ください。
欠損値を含む、先ほどの架空のtibbleを思い出してください。
R
rna_with_missing_values
出力
# A tibble: 7 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Apod GSM2545336 36194
3 Asl GSM2545337 361
4 Apod GSM2545337 10347
5 Asl GSM2545338 400
6 Apod GSM2545338 9173
7 Cyp2d22 GSM2545338 1603
R
wide_with_NA <- rna_with_missing_values |>
pivot_wider(names_from = sample,
values_from = expression)
wide_with_NA
出力
# A tibble: 3 × 4
gene GSM2545336 GSM2545337 GSM2545338
<chr> <dbl> <dbl> <dbl>
1 Asl 1170 361 400
2 Apod 36194 10347 9173
3 Cyp2d22 NA NA 1603
R
wide_with_NA |>
pivot_longer(names_to = "sample",
values_to = "expression",
-gene)
出力
# A tibble: 9 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Asl GSM2545337 361
3 Asl GSM2545338 400
4 Apod GSM2545336 36194
5 Apod GSM2545337 10347
6 Apod GSM2545338 9173
7 Cyp2d22 GSM2545336 NA
8 Cyp2d22 GSM2545337 NA
9 Cyp2d22 GSM2545338 1603
ワイド形式とロング形式へのピボットは、すべてのレプリケートが同じ構成になるようデータセットのバランスを整える上で有用な手法となります。
質問
rnaテーブルを起点として、pivot_wider()関数を用いて各マウスの遺伝子発現レベルを示すワイド形式のテーブルを作成してください。次に、pivot_longer()関数を用いてロング形式のテーブルに戻してください。
R
rna1 <- rna |>
select(gene, mouse, expression) |>
pivot_wider(names_from = mouse, values_from = expression)
rna1
出力
# A tibble: 1,474 × 23
gene `14` `9` `10` `15` `18` `6` `5` `11` `22` `13` `23`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Asl 1170 361 400 586 626 988 836 535 586 597 938
2 Apod 36194 10347 9173 10620 13021 29594 24959 13668 13230 15868 27769
3 Cyp2d22 4060 1616 1603 1901 2171 3349 3122 2008 2254 2277 2985
4 Klk6 287 629 641 578 448 195 186 1101 537 567 327
5 Fcrls 85 233 244 237 180 38 68 375 199 177 89
6 Slc2a4 782 231 248 265 313 786 528 249 266 357 654
7 Exd2 1619 2288 2235 2513 2366 1359 1474 3126 2379 2173 1531
8 Gjc2 288 595 568 551 310 146 186 791 454 370 240
9 Plp1 43217 101241 96534 58354 53126 27173 28728 98658 61356 61647 38019
10 Gnb4 1071 1791 1867 1430 1355 798 806 2437 1394 1554 960
# ℹ 1,464 more rows
# ℹ 11 more variables: `24` <dbl>, `8` <dbl>, `7` <dbl>, `1` <dbl>, `16` <dbl>,
# `21` <dbl>, `4` <dbl>, `2` <dbl>, `20` <dbl>, `12` <dbl>, `19` <dbl>
R
rna1 |>
pivot_longer(names_to = "mouse_id", values_to = "counts", -gene)
出力
# A tibble: 32,428 × 3
gene mouse_id counts
<chr> <chr> <dbl>
1 Asl 14 1170
2 Asl 9 361
3 Asl 10 400
4 Asl 15 586
5 Asl 18 626
6 Asl 6 988
7 Asl 5 836
8 Asl 11 535
9 Asl 22 586
10 Asl 13 597
# ℹ 32,418 more rows
質問
rnaデータフレームからX染色体およびY染色体に位置する遺伝子をサブセットし、sexを列、chromosome_nameを行、各染色体に位置する遺伝子の平均発現量を値として、以下のtibbleのようにデータフレームを展開します。

再整形する前に、集計する必要があります!
まずは、オスとメスのサンプルからX連鎖遺伝子およびY連鎖遺伝子の平均発現量を算出しましょう…
R
rna |>
filter(chromosome_name == "Y" | chromosome_name == "X") |>
group_by(sex, chromosome_name) |>
summarise(mean = mean(expression))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by sex and chromosome_name.
ℹ Output is grouped by sex.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(sex, chromosome_name))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 4 × 3
# Groups: sex [2]
sex chromosome_name mean
<chr> <chr> <dbl>
1 Female X 3504.
2 Female Y 3
3 Male X 2497.
4 Male Y 2117.
また、テーブルをワイド形式にピボットします。
R
rna_1 <- rna |>
filter(chromosome_name == "Y" | chromosome_name == "X") |>
group_by(sex, chromosome_name) |>
summarise(mean = mean(expression)) |>
pivot_wider(names_from = sex,
values_from = mean)
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by sex and chromosome_name.
ℹ Output is grouped by sex.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(sex, chromosome_name))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
R
rna_1
出力
# A tibble: 2 × 3
chromosome_name Female Male
<chr> <dbl> <dbl>
1 X 3504. 2497.
2 Y 3 2117.
次に、そのデータフレームをpivot_longer()で変換し、各行がchromosome_nameとgenderの一意の組み合わせとなるようにします。
R
rna_1 |>
pivot_longer(names_to = "gender",
values_to = "mean",
-chromosome_name)
出力
# A tibble: 4 × 3
chromosome_name gender mean
<chr> <chr> <dbl>
1 X Female 3504.
2 X Male 2497.
3 Y Female 3
4 Y Male 2117.
質問
rnaデータセットを使用して、各行が遺伝子の平均発現量を、各列が異なる時点を表す発現行列を作成してください。
まずは、遺伝子ごとおよび時間ごとの平均発現量を算出しましょう。
R
rna |>
group_by(gene, time) |>
summarise(mean_exp = mean(expression))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 4,422 × 3
# Groups: gene [1,474]
gene time mean_exp
<chr> <dbl> <dbl>
1 AI504432 0 1034.
2 AI504432 4 1104.
3 AI504432 8 1014
4 AW046200 0 155.
5 AW046200 4 152.
6 AW046200 8 81
7 AW551984 0 238
8 AW551984 4 302.
9 AW551984 8 342.
10 Aamp 0 4603.
# ℹ 4,412 more rows
pivot_wider()関数を使用する前に
R
rna_time <- rna |>
group_by(gene, time) |>
summarise(mean_exp = mean(expression)) |>
pivot_wider(names_from = time,
values_from = mean_exp)
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
R
rna_time
出力
# A tibble: 1,474 × 4
# Groups: gene [1,474]
gene `0` `4` `8`
<chr> <dbl> <dbl> <dbl>
1 AI504432 1034. 1104. 1014
2 AW046200 155. 152. 81
3 AW551984 238 302. 342.
4 Aamp 4603. 4870 4763.
5 Abca12 5.29 4.25 4.14
6 Abcc8 2576. 2609. 2292.
7 Abhd14a 591. 547. 432.
8 Abi2 4881. 4903. 4945.
9 Abi3bp 1175. 1061. 762.
10 Abl2 2170. 2078. 2131.
# ℹ 1,464 more rows
なお、これにより一部の列名が数字で始まるtibbleが生成されます。タイムポイントに対応する列を選択したい場合、列名をそのまま使うことはできません……。列4を選択するとどうなるでしょうか?
R
rna |>
group_by(gene, time) |>
summarise(mean_exp = mean(expression)) |>
pivot_wider(names_from = time,
values_from = mean_exp) |>
select(gene, 4)
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 1,474 × 2
# Groups: gene [1,474]
gene `8`
<chr> <dbl>
1 AI504432 1014
2 AW046200 81
3 AW551984 342.
4 Aamp 4763.
5 Abca12 4.14
6 Abcc8 2292.
7 Abhd14a 432.
8 Abi2 4945.
9 Abi3bp 762.
10 Abl2 2131.
# ℹ 1,464 more rows
タイムポイント4を選択するには、列名をバッククォート「`」で囲んで指定する必要があります。
R
rna |>
group_by(gene, time) |>
summarise(mean_exp = mean(expression)) |>
pivot_wider(names_from = time,
values_from = mean_exp) |>
select(gene, `4`)
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 1,474 × 2
# Groups: gene [1,474]
gene `4`
<chr> <dbl>
1 AI504432 1104.
2 AW046200 152.
3 AW551984 302.
4 Aamp 4870
5 Abca12 4.25
6 Abcc8 2609.
7 Abhd14a 547.
8 Abi2 4903.
9 Abi3bp 1061.
10 Abl2 2078.
# ℹ 1,464 more rows
別の方法としては、列名を数字で始まらない名前に変更することも考えられます。
R
rna |>
group_by(gene, time) |>
summarise(mean_exp = mean(expression)) |>
pivot_wider(names_from = time,
values_from = mean_exp) |>
rename("time0" = `0`, "time4" = `4`, "time8" = `8`) |>
select(gene, time4)
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
出力
# A tibble: 1,474 × 2
# Groups: gene [1,474]
gene time4
<chr> <dbl>
1 AI504432 1104.
2 AW046200 152.
3 AW551984 302.
4 Aamp 4870
5 Abca12 4.25
6 Abcc8 2609.
7 Abhd14a 547.
8 Abi2 4903.
9 Abi3bp 1061.
10 Abl2 2078.
# ℹ 1,464 more rows
質問
時点ごとの平均発現量を含む前回のデータフレームを使用し、時点8と時点0のフォールドチェンジ、および時点8と時点4のフォールドチェンジを格納する新しい列を作成してください。このテーブルを、算出されたフォールドチェンジをまとめたロングフォーマットのテーブルに変換してください。
rna_time ティブルから始めます:
R
rna_time
出力
# A tibble: 1,474 × 4
# Groups: gene [1,474]
gene `0` `4` `8`
<chr> <dbl> <dbl> <dbl>
1 AI504432 1034. 1104. 1014
2 AW046200 155. 152. 81
3 AW551984 238 302. 342.
4 Aamp 4603. 4870 4763.
5 Abca12 5.29 4.25 4.14
6 Abcc8 2576. 2609. 2292.
7 Abhd14a 591. 547. 432.
8 Abi2 4881. 4903. 4945.
9 Abi3bp 1175. 1061. 762.
10 Abl2 2170. 2078. 2131.
# ℹ 1,464 more rows
フォールドチェンジを算出します:
R
rna_time |>
mutate(time_8_vs_0 = `8` / `0`, time_8_vs_4 = `8` / `4`)
出力
# A tibble: 1,474 × 6
# Groups: gene [1,474]
gene `0` `4` `8` time_8_vs_0 time_8_vs_4
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 AI504432 1034. 1104. 1014 0.981 0.918
2 AW046200 155. 152. 81 0.522 0.532
3 AW551984 238 302. 342. 1.44 1.13
4 Aamp 4603. 4870 4763. 1.03 0.978
5 Abca12 5.29 4.25 4.14 0.784 0.975
6 Abcc8 2576. 2609. 2292. 0.889 0.878
7 Abhd14a 591. 547. 432. 0.731 0.791
8 Abi2 4881. 4903. 4945. 1.01 1.01
9 Abi3bp 1175. 1061. 762. 0.649 0.719
10 Abl2 2170. 2078. 2131. 0.982 1.03
# ℹ 1,464 more rows
また、pivot_longer() 関数を使用してください。
R
rna_time |>
mutate(time_8_vs_0 = `8` / `0`, time_8_vs_4 = `8` / `4`) |>
pivot_longer(names_to = "comparisons",
values_to = "Fold_changes",
time_8_vs_0:time_8_vs_4)
出力
# A tibble: 2,948 × 6
# Groups: gene [1,474]
gene `0` `4` `8` comparisons Fold_changes
<chr> <dbl> <dbl> <dbl> <chr> <dbl>
1 AI504432 1034. 1104. 1014 time_8_vs_0 0.981
2 AI504432 1034. 1104. 1014 time_8_vs_4 0.918
3 AW046200 155. 152. 81 time_8_vs_0 0.522
4 AW046200 155. 152. 81 time_8_vs_4 0.532
5 AW551984 238 302. 342. time_8_vs_0 1.44
6 AW551984 238 302. 342. time_8_vs_4 1.13
7 Aamp 4603. 4870 4763. time_8_vs_0 1.03
8 Aamp 4603. 4870 4763. time_8_vs_4 0.978
9 Abca12 5.29 4.25 4.14 time_8_vs_0 0.784
10 Abca12 5.29 4.25 4.14 time_8_vs_4 0.975
# ℹ 2,938 more rows
テーブルの結合
実際の場面では、データが複数のテーブルに分散していることがよくあります。通常、これは異なる種類の情報が異なる情報源から収集されるために起こります。
一部の解析では、すべてのテーブルに共通する列に基づいて、2つ以上のテーブルのデータを単一のデータフレームに統合することが望ましい場合があります。
dplyr
パッケージは、指定した列の一致に基づいて2つのデータフレームを結合するための一連の結合関数を提供します。ここでは、結合について簡単に解説します。詳細については、テーブル結合に関する章をご参照ください。データ変換チートシートでも、テーブル結合の概要を簡潔に説明しています。
元のrnaテーブルから3つの列と10行のみを抽出して作成する小さなテーブルrna_miniを用いて、結合の仕組みを説明します。
R
rna_mini <- rna |>
select(gene, sample, expression) |>
head(10)
rna_mini
出力
# A tibble: 10 × 3
gene sample expression
<chr> <chr> <dbl>
1 Asl GSM2545336 1170
2 Apod GSM2545336 36194
3 Cyp2d22 GSM2545336 4060
4 Klk6 GSM2545336 287
5 Fcrls GSM2545336 85
6 Slc2a4 GSM2545336 782
7 Exd2 GSM2545336 1619
8 Gjc2 GSM2545336 288
9 Plp1 GSM2545336 43217
10 Gnb4 GSM2545336 1071
2つ目のテーブル annot1
には、「gene」と「gene_description」の2つの列があります。リンクをクリックして
annot1.csv
をダウンロード し、data/
フォルダに移動させるか、あるいは以下のRコードを使用してフォルダに直接ダウンロードすることもできます。
R
download.file(url = "https://carpentries-incubator.github.io/bioc-intro/data/annot1.csv",
destfile = "data/annot1.csv")
annot1 <- read_csv(file = "data/annot1.csv")
annot1
出力
# A tibble: 10 × 2
gene gene_description
<chr> <chr>
1 Cyp2d22 cytochrome P450, family 2, subfamily d, polypeptide 22 [Source:MGI S…
2 Klk6 kallikrein related-peptidase 6 [Source:MGI Symbol;Acc:MGI:1343166]
3 Fcrls Fc receptor-like S, scavenger receptor [Source:MGI Symbol;Acc:MGI:19…
4 Plp1 proteolipid protein (myelin) 1 [Source:MGI Symbol;Acc:MGI:97623]
5 Exd2 exonuclease 3'-5' domain containing 2 [Source:MGI Symbol;Acc:MGI:192…
6 Apod apolipoprotein D [Source:MGI Symbol;Acc:MGI:88056]
7 Gnb4 guanine nucleotide binding protein (G protein), beta 4 [Source:MGI S…
8 Slc2a4 solute carrier family 2 (facilitated glucose transporter), member 4 …
9 Asl argininosuccinate lyase [Source:MGI Symbol;Acc:MGI:88084]
10 Gjc2 gap junction protein, gamma 2 [Source:MGI Symbol;Acc:MGI:2153060]
次に、dplyrパッケージのfull_join()関数を使用して、これら2つのテーブルをすべての変数を含む1つのテーブルに結合します。この関数は、1つ目のテーブルと2つ目のテーブルの列を照合するための共通変数を自動的に検出します。今回の例では、geneが共通変数となります。このような変数はキーと呼ばれます。キーは、異なるテーブル間で観測値を照合するために使用されます。
R
full_join(rna_mini, annot1)
出力
Joining with `by = join_by(gene)`
出力
# A tibble: 10 × 4
gene sample expression gene_description
<chr> <chr> <dbl> <chr>
1 Asl GSM2545336 1170 argininosuccinate lyase [Source:MGI Symbol;Acc…
2 Apod GSM2545336 36194 apolipoprotein D [Source:MGI Symbol;Acc:MGI:88…
3 Cyp2d22 GSM2545336 4060 cytochrome P450, family 2, subfamily d, polype…
4 Klk6 GSM2545336 287 kallikrein related-peptidase 6 [Source:MGI Sym…
5 Fcrls GSM2545336 85 Fc receptor-like S, scavenger receptor [Source…
6 Slc2a4 GSM2545336 782 solute carrier family 2 (facilitated glucose t…
7 Exd2 GSM2545336 1619 exonuclease 3'-5' domain containing 2 [Source:…
8 Gjc2 GSM2545336 288 gap junction protein, gamma 2 [Source:MGI Symb…
9 Plp1 GSM2545336 43217 proteolipid protein (myelin) 1 [Source:MGI Sym…
10 Gnb4 GSM2545336 1071 guanine nucleotide binding protein (G protein)…
実際には、遺伝子アノテーションの表記が異なる場合があります。
annot2
テーブルは、遺伝子名を含む変数のラベルが異なる点を除けば、annot1
と全く同じです。ここでも、ご自身で annot2.csv
をダウンロードして data/
に移動させるか、以下のRコードをご利用ください。
R
download.file(url = "https://carpentries-incubator.github.io/bioc-intro/data/annot2.csv",
destfile = "data/annot2.csv")
annot2 <- read_csv(file = "data/annot2.csv")
annot2
出力
# A tibble: 10 × 2
external_gene_name description
<chr> <chr>
1 Cyp2d22 cytochrome P450, family 2, subfamily d, polypeptide 22 [S…
2 Klk6 kallikrein related-peptidase 6 [Source:MGI Symbol;Acc:MGI…
3 Fcrls Fc receptor-like S, scavenger receptor [Source:MGI Symbol…
4 Plp1 proteolipid protein (myelin) 1 [Source:MGI Symbol;Acc:MGI…
5 Exd2 exonuclease 3'-5' domain containing 2 [Source:MGI Symbol;…
6 Apod apolipoprotein D [Source:MGI Symbol;Acc:MGI:88056]
7 Gnb4 guanine nucleotide binding protein (G protein), beta 4 [S…
8 Slc2a4 solute carrier family 2 (facilitated glucose transporter)…
9 Asl argininosuccinate lyase [Source:MGI Symbol;Acc:MGI:88084]
10 Gjc2 gap junction protein, gamma 2 [Source:MGI Symbol;Acc:MGI:…
変数名がいずれも一致しない場合は、照合に使用する変数を手動で指定できます。これらの変数は、以下に示すrna_miniおよびannot2テーブルのように、by引数を用いて指定できます。
R
full_join(rna_mini, annot2, by = c("gene" = "external_gene_name"))
出力
# A tibble: 10 × 4
gene sample expression description
<chr> <chr> <dbl> <chr>
1 Asl GSM2545336 1170 argininosuccinate lyase [Source:MGI Symbol;Acc…
2 Apod GSM2545336 36194 apolipoprotein D [Source:MGI Symbol;Acc:MGI:88…
3 Cyp2d22 GSM2545336 4060 cytochrome P450, family 2, subfamily d, polype…
4 Klk6 GSM2545336 287 kallikrein related-peptidase 6 [Source:MGI Sym…
5 Fcrls GSM2545336 85 Fc receptor-like S, scavenger receptor [Source…
6 Slc2a4 GSM2545336 782 solute carrier family 2 (facilitated glucose t…
7 Exd2 GSM2545336 1619 exonuclease 3'-5' domain containing 2 [Source:…
8 Gjc2 GSM2545336 288 gap junction protein, gamma 2 [Source:MGI Symb…
9 Plp1 GSM2545336 43217 proteolipid protein (myelin) 1 [Source:MGI Sym…
10 Gnb4 GSM2545336 1071 guanine nucleotide binding protein (G protein)…
上記からわかるように、結合後のテーブルでも最初のテーブルの変数名が保持されています。
課題:
こちらのリンクをクリックしてannot3テーブルをダウンロードし、data/リポジトリに保存してください。full_join()関数を使用して、テーブルrna_miniとannot3を結合してください。遺伝子Klk6、mt-Tf、mt-Rnr1、mt-Tv、mt-Rnr2、およびmt-Tl1にはどのような変化が見られましたか?
R
annot3 <- read_csv("data/annot3.csv")
full_join(rna_mini, annot3)
出力
# A tibble: 15 × 4
gene sample expression gene_description
<chr> <chr> <dbl> <chr>
1 Asl GSM2545336 1170 argininosuccinate lyase [Source:MGI Symbol;Acc…
2 Apod GSM2545336 36194 apolipoprotein D [Source:MGI Symbol;Acc:MGI:88…
3 Cyp2d22 GSM2545336 4060 cytochrome P450, family 2, subfamily d, polype…
4 Klk6 GSM2545336 287 <NA>
5 Fcrls GSM2545336 85 Fc receptor-like S, scavenger receptor [Source…
6 Slc2a4 GSM2545336 782 solute carrier family 2 (facilitated glucose t…
7 Exd2 GSM2545336 1619 exonuclease 3'-5' domain containing 2 [Source:…
8 Gjc2 GSM2545336 288 gap junction protein, gamma 2 [Source:MGI Symb…
9 Plp1 GSM2545336 43217 proteolipid protein (myelin) 1 [Source:MGI Sym…
10 Gnb4 GSM2545336 1071 guanine nucleotide binding protein (G protein)…
11 mt-Tf <NA> NA mitochondrially encoded tRNA phenylalanine [So…
12 mt-Rnr1 <NA> NA mitochondrially encoded 12S rRNA [Source:MGI S…
13 mt-Tv <NA> NA mitochondrially encoded tRNA valine [Source:MG…
14 mt-Rnr2 <NA> NA mitochondrially encoded 16S rRNA [Source:MGI S…
15 mt-Tl1 <NA> NA mitochondrially encoded tRNA leucine 1 [Source…
遺伝子 Klk6 は rna_mini にのみ存在し、遺伝子
mt-Tf、mt-Rnr1、mt-Tv、mt-Rnr2、および
mt-Tl1 は annot3
テーブルにのみ存在します。テーブルの各変数におけるこれらの値は、欠損値として符号化されています。
データのエクスポート
dplyrを用いて生データから情報を抽出したり要約したりする方法を習得したところで、これらの新しいデータセットをエクスポートして、共同研究者と共有したりアーカイブしたりしたい場合もあるでしょう。
CSVファイルをRに読み込む際に使用するread_csv()関数と同様に、データフレームからCSVファイルを生成するwrite_csv()関数があります。
write_csv()を使用する前に、作業ディレクトリ内に、生成されたデータセットを保存するための新しいフォルダdata_outputを作成します。生成されたデータセットを生データと同じディレクトリに保存することは避けるべきです。これらを分けて管理するのが良識です。dataフォルダには未加工の生データのみを格納し、削除や変更を行わないようそのままにしておきます。対照的に、スクリプトによってdata_outputディレクトリの内容が生成されるため、そこに含まれるファイルが削除されたとしても、いつでも再生成することが可能です。
write_csv()を使って、先ほど作成したrna_wideテーブルを保存しましょう。
R
write_csv(rna_wide, file = "data_output/rna_wide.csv")
- tidyverseメタパッケージを用いたRでの表形式データ
Content from データの可視化
最終更新日:2026-11-07 | ページの編集
所要時間: 120分
概要
質問
- Rによる可視化
目的
- ggplotを用いて、散布図、箱ひげ図、折れ線グラフなどを作成します。
- 共通のプロット設定を行います。
- ファセットとは何かを説明し、ggplotでファセットを適用してください。
- 既存のggplotプロットの外観(軸ラベルや色を含む)を変更してください。
- データフレーム内のデータから、複雑でカスタマイズされたプロットを作成します。
このエピソードは、Data Carpentriesの「生態学者のためのRによるデータ分析と可視化」レッスンに基づいています。
データの可視化
まず、必要なパッケージを読み込みます。ggplot2はtidyverseパッケージに含まれています。
R
library("tidyverse")
まだワークスペースにない場合は、前回のレッスンで保存したデータを読み込んでください。
R
rna <- read.csv("data/rnaseq.csv")
データ可視化チートシートでは、ggplot2の基本機能から高度な機能までを解説します。また、本パッケージで利用可能な多彩なデータ表現の概要を把握するとともに、備忘録としてもご活用いただけます。Thomas
Lin Pedersen氏による以下のビデオチュートリアル(パート1およびパート2)も非常に参考になります。
ggplot2によるプロット
ggplot2は、データフレーム内のデータから複雑なプロットを簡単に作成できるプロットパッケージです。どの変数をプロットするか、どのように表示するか、および全般的な視覚的プロパティを指定するための、よりプログラム的なインターフェースを提供します。ggplot2を支える理論的基盤は、Grammar
of Graphics( Wilkinson (2005)
)です。この手法を用いれば、基礎となるデータが変更された場合や、棒グラフから散布図に変更する場合でも、最小限の変更で済みます。これにより、最小限の調整や微調整で、論文掲載レベルのプロットを作成することが可能になります。
ggplot2 ( Wickham (2016) )
に関する書籍があり、概要を把握するのに適していますが、内容は古くなっています。第3版が準備中で、オンラインで無料で公開される予定です。ggplot2のウェブページ
(https://ggplot2.tidyverse.org)
には充実したドキュメントが掲載されています。
ggplot2
は、「ロング」形式(つまり、次元ごとに1列、観測値ごとに1行)のデータを想定して動作します。適切に構造化されたデータを使用することで、ggplot2
で図を作成する際の時間を大幅に短縮できます。
ggplotのグラフィックスは、新しい要素を順次追加することで構築されます。このようにレイヤーを重ねていくことで、プロットの自由度とカスタマイズ性を大幅に高めることができます。
Grammar of Graphicsの考え方は、あらゆるグラフを次の3つの構成要素から構築できるというものです。(1) データセット、(2) 座標系、(3) geoms — すなわち、データポイントを表す視覚的なマーク 1
ggplotを作成する際は、さまざまな種類のプロットに活用できる以下の基本テンプレートを使用します。
ggplot(data = <DATA>, mapping = aes(<MAPPINGS>)) + <GEOM_FUNCTION>()
-
ggplot()関数を使用し、data引数を用いてプロットを特定のデータフレームに紐付けます。
R
ggplot(data = rna)
- プロットする変数を選択し、それらをグラフ上でどのように表示するか(例:x/y座標、あるいはサイズ、形状、色などの属性)を指定することで、マッピング(aesthetic
(
aes) 関数を使用)を定義します。
R
ggplot(data = rna, mapping = aes(x = expression))
-
「geoms」――プロットにおけるデータの幾何学的表現(点、線、棒)を追加します。
ggplot2では多種多様なgeomsが提供されています。今回は、以下を含む代表的なものをいくつか使用します:* `geom_point()` for scatter plots, dot plots, etc. * `geom_histogram()` for histograms * `geom_boxplot()` for, well, boxplots! * `geom_line()` for trend lines, time series, etc.
プロットに図形を追加するには、+
演算子を使用します。まずは geom_histogram()
を使ってみましょう。
R
ggplot(data = rna, mapping = aes(x = expression)) +
geom_histogram()
出力
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

ggplot2パッケージの+は、既存のggplotオブジェクトを変更できるため、特に便利です。これにより、プロットテンプレートを簡単に作成し、さまざまな種類のプロットを手軽に探索できるようになります。したがって、上記のプロットは次のようなコードでも生成可能です。
R
# Assign plot to a variable
rna_plot <- ggplot(data = rna,
mapping = aes(x = expression))
# Draw the plot
rna_plot + geom_histogram()
課題
ヒストグラムを描画する際に、自動的にメッセージが表示されることにお気づきかと思います。
出力
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
geom_histogram() の引数 bins または
binwidth を変更することで、ビンの数や幅を変更できます。
R
# change bins
ggplot(rna, aes(x = expression)) +
geom_histogram(bins = 15)

R
# change binwidth
ggplot(rna, aes(x = expression)) +
geom_histogram(binwidth = 2000)

ここでは、データが右に歪んでいることが分かります。より対称的な分布にするために、底が2の対数変換を適用できます。なお、発現値が0の場合に-Infが返されるのを避けるため、ここでは小さな定数値(+1)を加算しています。
R
rna <- rna |>
mutate(expression_log = log2(expression + 1))
ここで、log2変換した発現量のヒストグラムを描くと、分布は確かに正規分布に近づいています。
R
ggplot(rna, aes(x = expression_log)) + geom_histogram()
出力
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

以降は、対数変換済みの発現値を扱います。
課題
この変換を可視化するもう一つの方法は、観測値のスケールを検討することです。例えば、プロット空間内で観測値をより適切に配置するために、軸のスケールを変更することが有効かもしれません。軸のスケール変更は、他のコンポーネントを追加・修正する場合と同様に(つまり、コマンドを段階的に追加することで)行います。実際にこの修正を試してみてください:
- 変換前の発現量をlog10スケールで表示してください。
scale_x_log10()を参照してください。先ほどのグラフと比較してください。なぜ今度は警告メッセージが表示されるようになったのでしょうか?
R
ggplot(data = rna,mapping = aes(x = expression))+
geom_histogram() +
scale_x_log10()
警告
Warning in scale_x_log10(): log-10 transformation introduced infinite values.
出力
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
警告
Warning: Removed 507 rows containing non-finite outside the scale range
(`stat_bin()`).

注記
-
ggplot()関数に記述した内容は、追加するあらゆるgeomレイヤーで反映されます(つまり、これらはプロット全体の設定となります)。これには、aes()で設定したx軸およびy軸のマッピングも含まれます。 - また、
ggplot()関数でグローバルに定義されたマッピングとは独立して、特定のgeomに対するマッピングを指定することも可能です。 - 新しいレイヤーを追加する際に使用する
+記号は、直前のレイヤーが含まれる行の末尾に配置する必要があります。代わりに、新しいレイヤーが含まれる行の先頭に+記号を付加した場合、ggplot2は新しいレイヤーを追加せず、エラーメッセージを返します。
R
# This is the correct syntax for adding layers
rna_plot +
geom_histogram()
# This will not add the new layer and will return an error message
rna_plot
+ geom_histogram()
プロットの反復的な作成
次に、2つの連続変数とgeom_point()関数を用いて散布図を作成します。このグラフでは、時刻8対時刻0、および時刻4対時刻0における発現量のlog2倍変化を表します。そのために、まず遺伝子ごと・時刻ごとの対数変換後の発現値の平均値を算出し、次に時刻8対時刻0、および時刻4対時刻0の平均対数発現値を差し引くことで、log倍変化を求めます。なお、ここでは後ほど遺伝子を表すために使用する遺伝子バイオタイプも含めています。倍変化の値は、rna_fc.という新しいデータフレームに保存します。
R
rna_fc <- rna |> select(gene, time,
gene_biotype, expression_log) |>
group_by(gene, time, gene_biotype) |>
summarize(mean_exp = mean(expression_log)) |>
pivot_wider(names_from = time,
values_from = mean_exp) |>
mutate(time_8_vs_0 = `8` - `0`, time_4_vs_0 = `4` - `0`)
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene, time, and gene_biotype.
ℹ Output is grouped by gene and time.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time, gene_biotype))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
次に、新たに作成したデータセット rna_fc を用いて ggplot
を作成します。ggplot2
を用いたプロット作成は、通常、反復的なプロセスとなります。まず、使用するデータセットを定義し、軸を配置し、geom
を選択します:
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
geom_point()

次に、このプロットを修正して、より多くの情報を抽出していきます。例えば、重なりを防ぐために透明度(alpha)を追加できます。
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
geom_point(alpha = 0.3)

また、すべての点に色を付けることもできます。
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
geom_point(alpha = 0.3, color = "blue")

あるいは、プロット内の各遺伝子を異なる色で塗り分けるには、引数
color
にベクトルを指定することもできます。ggplot2
は、ベクトル内の異なる値に対応するそれぞれ異なる色を割り当てます。以下に、gene_biotype
を用いて色分けを行う例を示します。
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
geom_point(alpha = 0.3, aes(color = gene_biotype))

また、ggplot()関数で指定されるマッピング内で直接色を指定することも可能です。この色はすべてのgeomレイヤーに反映され、マッピングはaes()で設定されたx軸およびy軸によって決定されます。
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0,
color = gene_biotype)) +
geom_point(alpha = 0.3)

最後に、geom_abline()関数を使って斜線を追加することもできます。
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0,
color = gene_biotype)) +
geom_point(alpha = 0.3) +
geom_abline(intercept = 0)

geomレイヤーをgeom_pointからgeom_jitterに変更しても、色は引き続きgene_biotypeによって決定される点にご留意ください。
R
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0,
color = gene_biotype)) +
geom_jitter(alpha = 0.3) +
geom_abline(intercept = 0)

課題
散布図は、小規模なデータセットにおいて有用な探索的ツールとなり得ます。rna_fcデータセットのように観測値数が多いデータセットでは、点の重なりが散布図の制約となる場合があります。このような状況に対処する一つの手法は、観測値を六角形ビンニングすることです。プロット領域は六角形に分割されます。各六角形には、その境界内に含まれる観測値の数に基づいて色が割り当てられます。
ggplot2で六角形ビンニングを使用するには、まずCRANからRパッケージhexbinをインストールし、ロードしてください。次に、
geom_hex()関数を使用してhexbin図を作成します。六角形ビンプロットは、散布図と比較してどのような相対的な長所と短所があるでしょうか。上記の散布図を確認し、作成した六角形ビンプロットと比較してください。
R
install.packages("hexbin")
R
library("hexbin")
ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
geom_hex() +
geom_abline(intercept = 0)

課題
先ほど学んだ内容を活用して、rnaデータセットからexpression_logをsampleに対してプロットした散布図を作成し、時間を異なる色で表示してください。このようなデータを示す方法として、これは適切でしょうか?
R
ggplot(data = rna, mapping = aes(y = expression_log, x = sample)) +
geom_point(aes(color = time))

箱ひげ図
箱ひげ図を用いて、各サンプルにおける遺伝子発現の分布を可視化できます。
R
ggplot(data = rna,
mapping = aes(y = expression_log, x = sample)) +
geom_boxplot()

箱ひげ図に点を追加することで、測定値の数やその分布をより明確に把握できます。
R
ggplot(data = rna,
mapping = aes(y = expression_log, x = sample)) +
geom_jitter(alpha = 0.2, color = "tomato") +
geom_boxplot(alpha = 0)

課題
ボックスプロットのレイヤーがジッターレイヤーよりも前面に配置されている点に注目してください。ボックスプロットを点の下に配置するには、コードのどこを変更すればよいでしょうか?
これら2つのgeomの順序を入れ替えるべきです。
R
ggplot(data = rna,
mapping = aes(y = expression_log, x = sample)) +
geom_boxplot(alpha = 0) +
geom_jitter(alpha = 0.2, color = "tomato")

x軸の値がまだ適切に読み取れないことに気づくかもしれません。ラベルの向きを変更し、重なり合わないように垂直方向と水平方向に調整しましょう。90度の角度を使うか、あるいは試行錯誤して斜め向きのラベルに適した角度を見つけてください:
R
ggplot(data = rna,
mapping = aes(y = expression_log, x = sample)) +
geom_jitter(alpha = 0.2, color = "tomato") +
geom_boxplot(alpha = 0) +
theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

課題
箱ひげ図のデータポイントに、感染期間(time)に応じた色を付けてください。
ヒント: クラス time を確認してください。ggplot
のマッピングで、time
のクラスを整数から因子に直接変更することを検討してください。なぜこれにより、R
のグラフ作成方法が変わるのでしょうか?
R
# time as integer
ggplot(data = rna,
mapping = aes(y = expression_log,
x = sample)) +
geom_jitter(alpha = 0.2, aes(color = time)) +
geom_boxplot(alpha = 0) +
theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

R
# time as factor
ggplot(data = rna,
mapping = aes(y = expression_log,
x = sample)) +
geom_jitter(alpha = 0.2, aes(color = as.factor(time))) +
geom_boxplot(alpha = 0) +
theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

課題
箱ひげ図は有用な要約手法ですが、分布の形状は見えなくなります。例えば、分布が二峰性である場合、箱ひげ図ではそれが分かりません。箱ひげ図に代わる手法としてバイオリンプロットがあり、こちらでは(点の密度の)形状が描かれます。
- 箱ひげ図をバイオリンプロットに置き換えてください。
geom_violin()を参照してください。引数fillを用いて、時間に応じてバイオリンプロットを塗りつぶしてください。
R
ggplot(data = rna,
mapping = aes(y = expression_log, x = sample)) +
geom_violin(aes(fill = as.factor(time))) +
theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

課題
- バイオリンプロットを修正し、
sexでバイオリンを塗りつぶしてください。
R
ggplot(data = rna,
mapping = aes(y = expression_log, x = sample)) +
geom_violin(aes(fill = sex)) +
theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

折れ線グラフ
まず、時刻0と時刻8を比較した際に対数フォールドチェンジが最大となった10個の遺伝子について、感染期間ごとの平均発現量を算出しましょう。まず、遺伝子を選択し、選択した10個の遺伝子を含むrnaのサブセットsub_rnaを作成します。次に、データをグループ化し、各グループ内での平均遺伝子発現量を算出します:
R
rna_fc <- rna_fc |> arrange(desc(time_8_vs_0))
genes_selected <- rna_fc$gene[1:10]
sub_rna <- rna |>
filter(gene %in% genes_selected)
mean_exp_by_time <- sub_rna |>
group_by(gene,time) |>
summarize(mean_exp = mean(expression_log))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene and time.
ℹ Output is grouped by gene.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
R
mean_exp_by_time
出力
# A tibble: 30 × 3
# Groups: gene [10]
gene time mean_exp
<chr> <int> <dbl>
1 Acr 0 5.07
2 Acr 4 5.54
3 Acr 8 7.31
4 Aipl1 0 3.70
5 Aipl1 4 3.89
6 Aipl1 8 6.56
7 Bst1 0 3.20
8 Bst1 4 3.77
9 Bst1 8 5.22
10 Chil3 0 4.00
# ℹ 20 more rows
x軸に感染期間、y軸に平均発現量をとった折れ線グラフを作成できます。
R
ggplot(data = mean_exp_by_time, mapping = aes(x = time, y = mean_exp)) +
geom_line()

残念ながら、すべての遺伝子のデータを一括してプロットしたため、これではうまくいきません。aesthetic関数にgroup = geneを含めるように修正し、ggplotに遺伝子ごとに線を描画させる必要があります。
R
ggplot(data = mean_exp_by_time,
mapping = aes(x = time, y = mean_exp, group = gene)) +
geom_line()

色を追加すれば、プロット上で遺伝子を区別できるようになります(color
を使用するとデータも自動的にグループ化されます):
R
ggplot(data = mean_exp_by_time,
mapping = aes(x = time, y = mean_exp, color = gene)) +
geom_line()

ファセット化
ggplot2には、データセットに含まれる因子に基づいて1つのプロットを複数の(サブ)プロットに分割できる「ファセット化」という特殊な手法があります。これらの異なるサブプロットは、直接比較しやすくするため、同一のプロパティ(軸の範囲、目盛りなど)を継承します。これを用いて、各遺伝子の経時変化を示す折れ線グラフを作成します。
R
ggplot(data = mean_exp_by_time,
mapping = aes(x = time, y = mean_exp)) + geom_line() +
facet_wrap(~ gene)

ここでは、すべてのサブプロットにおいてx軸とy軸のスケールが同一になっています。scalesを修正してy軸のスケールを自由に設定することで、このデフォルトの動作を変更できます。
R
ggplot(data = mean_exp_by_time,
mapping = aes(x = time, y = mean_exp)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y")

次に、各プロットの線をマウスの性別ごとに分けたいと思います。そのためには、データフレーム内でgene、time、およびsexでグループ化した平均発現量を算出する必要があります:
R
mean_exp_by_time_sex <- sub_rna |>
group_by(gene, time, sex) |>
summarize(mean_exp = mean(expression_log))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by gene, time, and sex.
ℹ Output is grouped by gene and time.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(gene, time, sex))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
R
mean_exp_by_time_sex
出力
# A tibble: 60 × 4
# Groups: gene, time [30]
gene time sex mean_exp
<chr> <int> <chr> <dbl>
1 Acr 0 Female 5.13
2 Acr 0 Male 5.00
3 Acr 4 Female 5.93
4 Acr 4 Male 5.15
5 Acr 8 Female 7.27
6 Acr 8 Male 7.36
7 Aipl1 0 Female 3.67
8 Aipl1 0 Male 3.73
9 Aipl1 4 Female 4.07
10 Aipl1 4 Male 3.72
# ℹ 50 more rows
次に、colorを用いて(単一のプロット内で)さらに性別で分割することで、面分けプロットを作成できます。
R
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y")

通常、背景が白のプロットは印刷した際により見やすくなります。theme_bw()関数を使用して、背景を白に設定できます。さらに、グリッドを削除することも可能です。
R
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank())

課題
先ほど学んだ内容を活用して、感染期間を通じて各染色体の平均発現量がどのように変化するかを示すプロットを作成してください。
R
mean_exp_by_chromosome <- rna |>
group_by(chromosome_name, time) |>
summarize(mean_exp = mean(expression_log))
出力
`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by chromosome_name and time.
ℹ Output is grouped by chromosome_name.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(chromosome_name, time))` for per-operation grouping
(`?dplyr::dplyr_by`) instead.
R
ggplot(data = mean_exp_by_chromosome, mapping = aes(x = time,
y = mean_exp)) +
geom_line() +
facet_wrap(~ chromosome_name, scales = "free_y")

facet_wrap
ジオメトリでは、プロットを任意の次元数で抽出できるため、1ページにすっきりと収めることができます。一方、facet_grid
ジオメトリでは、数式表記を用いてプロットの配置方法を明示的に指定できます(rows ~ columns;.
は、1行または1列のみを示すプレースホルダーとして使用できます)。
先ほどのプロットを修正して、オスとメスの平均遺伝子発現が経時的にどのように変化したかを比較してみましょう。
R
# One column, facet by rows
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = gene)) +
geom_line() +
facet_grid(sex ~ .)

R
# One row, facet by column
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = gene)) +
geom_line() +
facet_grid(. ~ sex)

ggplot2 テーマ
プロットの背景を白に変更するtheme_bw()のほか、ggplot2には可視化の外観を素早く変更する際に役立つテーマが他にもいくつか用意されています。テーマの全一覧はhttps://ggplot2.tidyverse.org/reference/ggtheme.htmlで確認できます。theme_minimal()とtheme_light()は人気があり、theme_void()は独自のテーマを新規に作成する際の出発点として有用です。
ggthemes
パッケージでは、多彩なオプション(Excel 2003
テーマを含む)が提供されています。ggplot2
拡張機能ウェブサイト では、追加テーマを含む、ggplot2
の機能を拡張するパッケージの一覧が掲載されています。
カスタマイズ
時間と遺伝子ごとの平均発現量を、性別で色分けしたファセットプロットに戻りましょう。
ggplot2
チートシートを確認し、プロットを改善する方法を考えてみてください。
次に、軸の名前を「time」や「mean_exp」よりも分かりやすいものに変更し、図にタイトルを追加します。
R
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank()) +
labs(title = "Mean gene expression by duration of the infection",
x = "Duration of the infection (in days)",
y = "Mean gene expression")

軸にはより分かりやすい名称が付けられていますが、フォントサイズを大きくすることで視認性を向上させることができます。
R
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank()) +
labs(title = "Mean gene expression by duration of the infection",
x = "Duration of the infection (in days)",
y = "Mean gene expression") +
theme(text = element_text(size = 16))

なお、プロットのフォントを変更することも可能です。Windowsをご利用の場合は、extrafontパッケージをインストールする必要がある場合があります。
さらに、x軸・y軸のテキストの色やグリッドの色などをカスタマイズできます。また、例えばlegend.positionを"top"に設定することで、凡例を上部に移動させることも可能です。
R
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank()) +
labs(title = "Mean gene expression by duration of the infection",
x = "Duration of the infection (in days)",
y = "Mean gene expression") +
theme(text = element_text(size = 16),
axis.text.x = element_text(colour = "royalblue4", size = 12),
axis.text.y = element_text(colour = "royalblue4", size = 12),
panel.grid = element_line(colour="lightsteelblue1"),
legend.position = "top")

作成した変更がデフォルトのテーマよりも気に入った場合は、オブジェクトとして保存しておけば、今後作成する他のプロットにも簡単に適用できます。以下に、以前に作成したヒストグラムを用いた例を示します。
R
blue_theme <- theme(axis.text.x = element_text(colour = "royalblue4",
size = 12),
axis.text.y = element_text(colour = "royalblue4",
size = 12),
text = element_text(size = 16),
panel.grid = element_line(colour="lightsteelblue1"))
ggplot(rna, aes(x = expression_log)) +
geom_histogram(bins = 20) +
blue_theme

課題
以上の情報を踏まえ、さらに5分ほど時間をかけて、この演習で生成したプロットのいずれかを改善するか、あるいはご自身で見栄えの良いグラフを作成してみてください。参考にRStudioのggplot2チートシートをご活用ください。いくつかアイデアを挙げます:
- 線の太さを変えられるか試してみてください。
- 凡例の名前を変更する方法はありますか?ラベルはどうでしょうか?(ヒント:
scale_で始まるggplot関数を探してみてください) - 別のカラーパレットを試すか、線の色を手動で指定してみてください(https://www.cookbook-r.com/Graphs/Colors_(ggplot2)/を参照してください)。
例えば、このプロットに基づくと:
R
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank())

以下のようにカスタマイズできます。
R
# change the thickness of the lines
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line(size=1.5) +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank())
警告
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once per session.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

R
# change the name of the legend and the labels
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank()) +
scale_color_discrete(name = "Gender", labels = c("F", "M"))

R
# using a different color palette
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank()) +
scale_color_brewer(name = "Gender", labels = c("F", "M"), palette = "Dark2")

R
# manually specifying the colors
ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
theme_bw() +
theme(panel.grid = element_blank()) +
scale_color_manual(name = "Gender", labels = c("F", "M"),
values = c("royalblue", "deeppink"))

プロットの合成
ファセット機能は、1つのプロットを複数のサブプロットに分割するのに便利なツールです。しかし、異なる変数や、さらには異なるデータフレームに基づいた複数の独立したプロットを、1つの図にまとめたい場合もあるでしょう。
まずは、並べて表示したい2つのプロットを作成しましょう。
最初のグラフでは、染色体ごとのユニークな遺伝子数をカウントしています。まず、chromosome_nameの水準を並べ替え、染色体ごとのユニークな遺伝子をフィルタリングします。また、視認性を高めるため、y軸のスケールをlog10スケールに変更します。
R
rna$chromosome_name <- factor(rna$chromosome_name,
levels = c(1:19,"X","Y"))
count_gene_chromosome <- rna |>
select(chromosome_name, gene) |>
distinct() |> ggplot() +
geom_bar(aes(x = chromosome_name), fill = "seagreen",
position = "dodge", stat = "count") +
labs(y = "log10(n genes)", x = "chromosome") +
scale_y_log10()
count_gene_chromosome

以下では、legend.positionを"none"に設定することで、凡例を完全に削除します。
R
exp_boxplot_sex <- ggplot(rna, aes(y=expression_log, x = as.factor(time),
color=sex)) +
geom_boxplot(alpha = 0) +
labs(y = "Mean gene exp",
x = "time") + theme(legend.position = "none")
exp_boxplot_sex

patchwork
パッケージは、+
を用いて図を(通常は左右に)配置することで、図を組み合わせるための洗練された手法を提供します。より具体的には、|
は図を明示的に左右に配置し、/
は図を上下に積み重ねます。
R
install.packages("patchwork")
R
library("patchwork")
count_gene_chromosome + exp_boxplot_sex

R
## or count_gene_chromosome | exp_boxplot_sex
R
count_gene_chromosome / exp_boxplot_sex

plot_layoutを用いて最終的な構成のレイアウトをさらに調整することで、より複雑なレイアウトを作成できます。
R
count_gene_chromosome + exp_boxplot_sex + plot_layout(ncol = 1)

R
count_gene_chromosome +
(count_gene_chromosome + exp_boxplot_sex) +
exp_boxplot_sex +
plot_layout(ncol = 1)

最後のプロットは、| と /
のコンポーザーを使用して作成することもできます。
R
count_gene_chromosome /
(count_gene_chromosome | exp_boxplot_sex) /
exp_boxplot_sex

patchworkの詳細については、ウェブページまたはこちらの動画をご覧ください。
別の選択肢として、grid.arrange()を用いて個別のggplotを1つの図に統合できるgridExtraパッケージがあります。
R
install.packages("gridExtra")
R
library("gridExtra")
grid.arrange(count_gene_chromosome, exp_boxplot_sex, ncol = 2)

シンプルな配置に用いるncolやnrowといった引数に加え、より複雑なレイアウトを構築するためのツールも用意されています。
プロットのエクスポート
プロットを作成した後、お好みの形式でファイルに保存できます。RStudioのPlotペインにある「Export」タブでは、プロットが低解像度で保存されます。この解像度は多くの学術誌で受理されず、ポスター用としても適切に拡大できません。
代わりに、ggsave()関数を使用してください。この関数では、適切な引数(width、height、dpi)を調整することで、プロットの寸法や解像度を簡単に変更できます。
作業ディレクトリにfig_output/フォルダがあることを確認してください。
R
my_plot <- ggplot(data = mean_exp_by_time_sex,
mapping = aes(x = time, y = mean_exp, color = sex)) +
geom_line() +
facet_wrap(~ gene, scales = "free_y") +
labs(title = "Mean gene expression by duration of the infection",
x = "Duration of the infection (in days)",
y = "Mean gene expression") +
guides(color=guide_legend(title="Gender")) +
theme_bw() +
theme(axis.text.x = element_text(colour = "royalblue4", size = 12),
axis.text.y = element_text(colour = "royalblue4", size = 12),
text = element_text(size = 16),
panel.grid = element_line(colour="lightsteelblue1"),
legend.position = "top")
ggsave("fig_output/mean_exp_by_time_sex.png", my_plot, width = 15,
height = 10)
# This also works for grid.arrange() plots
combo_plot <- grid.arrange(count_gene_chromosome, exp_boxplot_sex,
ncol = 2, widths = c(4, 6))
ggsave("fig_output/combo_plot_chromosome_sex.png", combo_plot,
width = 10, dpi = 300)
注:パラメータ width と height
は、保存されるプロットのフォントサイズも決定します。
可視化用のその他のパッケージ
ggplot2は非常に強力なパッケージであり、私たちの整然データおよび整然ツールのパイプラインに非常にうまく適合します。Rには、他にも見逃せない可視化パッケージがあります。
基本グラフィックス
Rに標準で搭載されているグラフィックスシステム(一般にbase R
graphicsと呼ばれます)は、シンプルで高速です。これはペインターモデルまたはキャンバスモデルに基づいており、異なる出力が直接重ね合わせられます(下の図参照)。この点は、専用のオブジェクトを返し、それらを画面やファイルにレンダリングし、さらには更新することさえ可能なggplot2(および後述するlattice)とは根本的に異なります。
R
par(mfrow = c(1, 3))
plot(1:20, main = "First layer, produced with plot(1:20)")
plot(1:20, main = "A horizontal red line, added with abline(h = 10)")
abline(h = 10, col = "red")
plot(1:20, main = "A rectangle, added with rect(5, 5, 15, 15)")
abline(h = 10, col = "red")
rect(5, 5, 15, 15, lwd = 3)

もう一つの主な違いは、base
graphicsのプロット関数が入力型に応じて適切な処理を行おうとする点です。つまり、入力のクラスに基づいて動作を適応させます。これは、データフレームのみを入力として受け付け、プロットを一つひとつ構築する必要があるggplot2の仕組みとは、ここでも大きく異なります。
R
par(mfrow = c(2, 2))
boxplot(rnorm(100),
main = "Boxplot of rnorm(100)")
boxplot(matrix(rnorm(100), ncol = 10),
main = "Boxplot of matrix(rnorm(100), ncol = 10)")
hist(rnorm(100))
hist(matrix(rnorm(100), ncol = 10))

baseグラフィックスの標準機能による手法は、シンプルで標準的な図表であれば非常に効率的であり、plotやhist、boxplotといった単一の関数とわずか1行のコードで極めて迅速に作成できます。しかし、デフォルト設定が必ずしも最適とは限らず、特に図表が複雑になる場合(例えばファセットを作成する場合など)、調整作業に多大な時間と手間がかかることがあります。
latticeパッケージ
lattice
パッケージは、データフレームを入力として用い、グラフィカルなオブジェクトを返し、ファセット機能をサポートするという点で
ggplot2 と似ています。一方、lattice
はグラフィックスの文法に基づいておらず、インターフェースもより複雑です。
latticeパッケージの優れたリファレンスは Deepayan (2008). です。
- Rによる可視化
出典:データ可視化チートシート。↩︎
Content from 次のステップ
最終更新日:2026-11-07 | ページの編集
所要時間: 90分
概要
質問
-
SummarizedExperimentとは何ですか? - Bioconductorとは何ですか?
目的
- Bioconductorプロジェクトについて紹介します。
- データコンテナの概念を紹介します。
- オミクス解析で広く用いられる
SummarizedExperimentの概要を説明してください。
次のステップ
バイオインフォマティクスにおけるデータは、複雑な場合が多くあります。これに対処するため、開発者は扱う必要のあるデータの特性に合わせた専用のデータコンテナ(クラスと呼ばれます)を定義します。
この側面は、パッケージ間で共通のコアデータ基盤を採用しているBioconductor1プロジェクトの中核をなすものです。これは間違いなくBioconductorの成功に寄与しました。Bioconductorのパッケージ開発者は、プロジェクト全体に一貫性、相互運用性、および安定性をもたらすため、既存の基盤を活用することが推奨されています。
このようなオミクスデータコンテナの例として、SummarizedExperimentクラスを紹介します。
SummarizedExperiment
下図は、SummarizedExperimentクラスの構成を示しています。
SummarizedExperimentクラスのオブジェクトには、以下が含まれます。
行列形式のオブジェクトとして格納された定量的なオミクスデータ(発現データ)を含む1つ(または複数)のアッセイです。行方向に特徴量(遺伝子、転写産物、タンパク質など)が、列方向にサンプルが定義されています。
サンプルの共変量を格納したサンプルメタデータスロットです。データフレームとして保存されています。このテーブルの行はサンプルを表します(行は発現データの列と完全に一致します)。
特徴量の共変量を格納した特徴量メタデータスロットです。データフレームとして保存されています。このデータフレームの行は、発現データの行と完全に一致します。
SummarizedExperimentの整合性により、データ操作時に異なるスロットの次元が常に一致することが保証されます(すなわち、発現データの列とサンプルメタデータの行、および発現データの行と特徴量メタデータの行が一致します)。例えば、アッセイから1つのサンプルを除外する場合、同一の操作によってサンプルメタデータからも自動的に削除されます。
メタデータスロットには、他の構造に影響を与えることなく、追加の共変量(列)を加えることができます。
SummarizedExperimentの作成
SummarizedExperimentを作成するため、個々の構成要素、すなわちカウント行列、サンプルおよび遺伝子のメタデータをCSVファイルから作成します。これらは通常、RNA-Seqデータ(生データの処理後)の提供形式となっています。
-
発現行列:カウント行列を読み込み、最初の列に行名/遺伝子名が含まれていることを指定して、
data.frameをmatrixに変換します。こちらのリンクからダウンロードいただけます。
R
count_matrix <- read.csv("data/count_matrix.csv",
row.names = 1) |>
as.matrix()
count_matrix[1:5, ]
出力
GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
Asl 1170 361 400 586 626 988
Apod 36194 10347 9173 10620 13021 29594
Cyp2d22 4060 1616 1603 1901 2171 3349
Klk6 287 629 641 578 448 195
Fcrls 85 233 244 237 180 38
GSM2545342 GSM2545343 GSM2545344 GSM2545345 GSM2545346 GSM2545347
Asl 836 535 586 597 938 1035
Apod 24959 13668 13230 15868 27769 34301
Cyp2d22 3122 2008 2254 2277 2985 3452
Klk6 186 1101 537 567 327 233
Fcrls 68 375 199 177 89 67
GSM2545348 GSM2545349 GSM2545350 GSM2545351 GSM2545352 GSM2545353
Asl 494 481 666 937 803 541
Apod 11258 11812 15816 29242 20415 13682
Cyp2d22 1883 2014 2417 3678 2920 2216
Klk6 742 881 828 250 798 710
Fcrls 300 233 231 81 303 285
GSM2545354 GSM2545362 GSM2545363 GSM2545380
Asl 473 748 576 1192
Apod 11088 15916 11166 38148
Cyp2d22 1821 2842 2011 4019
Klk6 894 501 598 259
Fcrls 248 179 184 68
R
dim(count_matrix)
出力
[1] 1474 22
- サンプルの説明表は、こちらのリンクからご覧いただけます。
R
sample_metadata <- read.csv("data/sample_metadata.csv")
sample_metadata
出力
sample organism age sex infection strain time tissue mouse
1 GSM2545336 Mus musculus 8 Female InfluenzaA C57BL/6 8 Cerebellum 14
2 GSM2545337 Mus musculus 8 Female NonInfected C57BL/6 0 Cerebellum 9
3 GSM2545338 Mus musculus 8 Female NonInfected C57BL/6 0 Cerebellum 10
4 GSM2545339 Mus musculus 8 Female InfluenzaA C57BL/6 4 Cerebellum 15
5 GSM2545340 Mus musculus 8 Male InfluenzaA C57BL/6 4 Cerebellum 18
6 GSM2545341 Mus musculus 8 Male InfluenzaA C57BL/6 8 Cerebellum 6
7 GSM2545342 Mus musculus 8 Female InfluenzaA C57BL/6 8 Cerebellum 5
8 GSM2545343 Mus musculus 8 Male NonInfected C57BL/6 0 Cerebellum 11
9 GSM2545344 Mus musculus 8 Female InfluenzaA C57BL/6 4 Cerebellum 22
10 GSM2545345 Mus musculus 8 Male InfluenzaA C57BL/6 4 Cerebellum 13
11 GSM2545346 Mus musculus 8 Male InfluenzaA C57BL/6 8 Cerebellum 23
12 GSM2545347 Mus musculus 8 Male InfluenzaA C57BL/6 8 Cerebellum 24
13 GSM2545348 Mus musculus 8 Female NonInfected C57BL/6 0 Cerebellum 8
14 GSM2545349 Mus musculus 8 Male NonInfected C57BL/6 0 Cerebellum 7
15 GSM2545350 Mus musculus 8 Male InfluenzaA C57BL/6 4 Cerebellum 1
16 GSM2545351 Mus musculus 8 Female InfluenzaA C57BL/6 8 Cerebellum 16
17 GSM2545352 Mus musculus 8 Female InfluenzaA C57BL/6 4 Cerebellum 21
18 GSM2545353 Mus musculus 8 Female NonInfected C57BL/6 0 Cerebellum 4
19 GSM2545354 Mus musculus 8 Male NonInfected C57BL/6 0 Cerebellum 2
20 GSM2545362 Mus musculus 8 Female InfluenzaA C57BL/6 4 Cerebellum 20
21 GSM2545363 Mus musculus 8 Male InfluenzaA C57BL/6 4 Cerebellum 12
22 GSM2545380 Mus musculus 8 Female InfluenzaA C57BL/6 8 Cerebellum 19
R
dim(sample_metadata)
出力
[1] 22 9
- 遺伝子に関する表は、こちらのリンクからご覧いただけます。
R
gene_metadata <- read.csv("data/gene_metadata.csv")
gene_metadata[1:10, 1:4]
出力
gene ENTREZID
1 Asl 109900
2 Apod 11815
3 Cyp2d22 56448
4 Klk6 19144
5 Fcrls 80891
6 Slc2a4 20528
7 Exd2 97827
8 Gjc2 118454
9 Plp1 18823
10 Gnb4 14696
product
1 argininosuccinate lyase, transcript variant X1
2 apolipoprotein D, transcript variant 3
3 cytochrome P450, family 2, subfamily d, polypeptide 22, transcript variant 2
4 kallikrein related-peptidase 6, transcript variant 2
5 Fc receptor-like S, scavenger receptor, transcript variant X1
6 solute carrier family 2 (facilitated glucose transporter), member 4
7 exonuclease 3'-5' domain containing 2
8 gap junction protein, gamma 2, transcript variant 1
9 proteolipid protein (myelin) 1, transcript variant 1
10 guanine nucleotide binding protein (G protein), beta 4, transcript variant X2
ensembl_gene_id
1 ENSMUSG00000025533
2 ENSMUSG00000022548
3 ENSMUSG00000061740
4 ENSMUSG00000050063
5 ENSMUSG00000015852
6 ENSMUSG00000018566
7 ENSMUSG00000032705
8 ENSMUSG00000043448
9 ENSMUSG00000031425
10 ENSMUSG00000027669
R
dim(gene_metadata)
出力
[1] 1474 9
これらのテーブルからSummarizedExperimentを作成します。
assayとして使用されるカウント行列サンプルを記述した表は、サンプルメタデータスロットとして使用されます。
遺伝子を記述したテーブルは、features metadataスロットとして使用されます。
これを行うには、SummarizedExperiment
コンストラクタを用いて各パーツを組み合わせることができます。
R
## BiocManager::install("SummarizedExperiment")
library("SummarizedExperiment")
まず、カウント行列とサンプルアノテーションにおいてサンプルの順序が一致していること、およびカウント行列と遺伝子アノテーションにおいて遺伝子の順序が一致していることを確認します。
R
stopifnot(rownames(count_matrix) == gene_metadata$gene)
stopifnot(colnames(count_matrix) == sample_metadata$sample)
R
se <- SummarizedExperiment(assays = list(counts = count_matrix),
colData = sample_metadata,
rowData = gene_metadata)
se
出力
class: SummarizedExperiment
dim: 1474 22
metadata(0):
assays(1): counts
rownames(1474): Asl Apod ... Lmx1a Pbx1
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(22): GSM2545336 GSM2545337 ... GSM2545363 GSM2545380
colData names(9): sample organism ... tissue mouse
データの保存
前回のエピソードで行ったようなスプレッドシートへのデータエクスポートには、第1章で述べた点(小数点区切り文字に関する,と.の不整合の可能性や、変数型の定義がないことなど)など、いくつかの制約があります。さらに、スプレッドシートへのデータエクスポートが有効なのは、データフレームや行列のような矩形データに限られます。
R特有の、いかなるオペレーティングシステムでも確実に動作するより汎用的なデータ保存方法は、saveRDS関数を使用することです。このようにオブジェクトを保存すると、ディスク上にバイナリ形式のデータ(ここではrdsというファイル拡張子を使用)が生成され、これをreadRDS関数を用いてRに読み込むことができます。
R
saveRDS(se, file = "data_output/se.rds")
rm(se)
se <- readRDS("data_output/se.rds")
head(se)
まとめとして、Rで保存し、再度Rで読み込むデータについては、saveRDSおよびreadRDSを用いて保存・読み込みを行う方法が推奨されます。表形式のデータをRを使用していない相手と共有する必要がある場合は、テキストベースのスプレッドシートにエクスポートすることも有効な選択肢となります。
このデータ構造を用いることで、assay関数で発現行列にアクセスできます。
R
head(assay(se))
出力
GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
Asl 1170 361 400 586 626 988
Apod 36194 10347 9173 10620 13021 29594
Cyp2d22 4060 1616 1603 1901 2171 3349
Klk6 287 629 641 578 448 195
Fcrls 85 233 244 237 180 38
Slc2a4 782 231 248 265 313 786
GSM2545342 GSM2545343 GSM2545344 GSM2545345 GSM2545346 GSM2545347
Asl 836 535 586 597 938 1035
Apod 24959 13668 13230 15868 27769 34301
Cyp2d22 3122 2008 2254 2277 2985 3452
Klk6 186 1101 537 567 327 233
Fcrls 68 375 199 177 89 67
Slc2a4 528 249 266 357 654 693
GSM2545348 GSM2545349 GSM2545350 GSM2545351 GSM2545352 GSM2545353
Asl 494 481 666 937 803 541
Apod 11258 11812 15816 29242 20415 13682
Cyp2d22 1883 2014 2417 3678 2920 2216
Klk6 742 881 828 250 798 710
Fcrls 300 233 231 81 303 285
Slc2a4 271 304 349 715 513 320
GSM2545354 GSM2545362 GSM2545363 GSM2545380
Asl 473 748 576 1192
Apod 11088 15916 11166 38148
Cyp2d22 1821 2842 2011 4019
Klk6 894 501 598 259
Fcrls 248 179 184 68
Slc2a4 248 350 317 796
R
dim(assay(se))
出力
[1] 1474 22
colData関数を使用して、サンプルのメタデータにアクセスできます。
R
colData(se)
出力
DataFrame with 22 rows and 9 columns
sample organism age sex infection
<character> <character> <integer> <character> <character>
GSM2545336 GSM2545336 Mus musculus 8 Female InfluenzaA
GSM2545337 GSM2545337 Mus musculus 8 Female NonInfected
GSM2545338 GSM2545338 Mus musculus 8 Female NonInfected
GSM2545339 GSM2545339 Mus musculus 8 Female InfluenzaA
GSM2545340 GSM2545340 Mus musculus 8 Male InfluenzaA
... ... ... ... ... ...
GSM2545353 GSM2545353 Mus musculus 8 Female NonInfected
GSM2545354 GSM2545354 Mus musculus 8 Male NonInfected
GSM2545362 GSM2545362 Mus musculus 8 Female InfluenzaA
GSM2545363 GSM2545363 Mus musculus 8 Male InfluenzaA
GSM2545380 GSM2545380 Mus musculus 8 Female InfluenzaA
strain time tissue mouse
<character> <integer> <character> <integer>
GSM2545336 C57BL/6 8 Cerebellum 14
GSM2545337 C57BL/6 0 Cerebellum 9
GSM2545338 C57BL/6 0 Cerebellum 10
GSM2545339 C57BL/6 4 Cerebellum 15
GSM2545340 C57BL/6 4 Cerebellum 18
... ... ... ... ...
GSM2545353 C57BL/6 0 Cerebellum 4
GSM2545354 C57BL/6 0 Cerebellum 2
GSM2545362 C57BL/6 4 Cerebellum 20
GSM2545363 C57BL/6 4 Cerebellum 12
GSM2545380 C57BL/6 8 Cerebellum 19
R
dim(colData(se))
出力
[1] 22 9
また、rowData関数を使用してフィーチャーのメタデータにアクセスすることもできます。
R
head(rowData(se))
出力
DataFrame with 6 rows and 9 columns
gene ENTREZID product ensembl_gene_id
<character> <integer> <character> <character>
Asl Asl 109900 argininosuccinate ly.. ENSMUSG00000025533
Apod Apod 11815 apolipoprotein D, tr.. ENSMUSG00000022548
Cyp2d22 Cyp2d22 56448 cytochrome P450, fam.. ENSMUSG00000061740
Klk6 Klk6 19144 kallikrein related-p.. ENSMUSG00000050063
Fcrls Fcrls 80891 Fc receptor-like S, .. ENSMUSG00000015852
Slc2a4 Slc2a4 20528 solute carrier famil.. ENSMUSG00000018566
external_synonym chromosome_name gene_biotype phenotype_description
<character> <character> <character> <character>
Asl 2510006M18Rik 5 protein_coding abnormal circulating..
Apod NA 16 protein_coding abnormal lipid homeo..
Cyp2d22 2D22 15 protein_coding abnormal skin morpho..
Klk6 Bssp 7 protein_coding abnormal cytokine le..
Fcrls 2810439C17Rik 3 protein_coding decreased CD8-positi..
Slc2a4 Glut-4 11 protein_coding abnormal circulating..
hsapiens_homolog_associated_gene_name
<character>
Asl ASL
Apod APOD
Cyp2d22 CYP2D6
Klk6 KLK6
Fcrls FCRL2
Slc2a4 SLC2A4
R
dim(rowData(se))
出力
[1] 1474 9
SummarizedExperimentのサブセット抽出
SummarizedExperimentは、データフレームと同様に、数値、文字、または論理値を用いてサブセットを抽出できます。
以下では、最初の3つのサンプルの最初の5つの特徴量のみを含むSummarizedExperimentクラスの新しいインスタンスを作成します。
R
se1 <- se[1:5, 1:3]
se1
出力
class: SummarizedExperiment
dim: 5 3
metadata(0):
assays(1): counts
rownames(5): Asl Apod Cyp2d22 Klk6 Fcrls
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(3): GSM2545336 GSM2545337 GSM2545338
colData names(9): sample organism ... tissue mouse
R
colData(se1)
出力
DataFrame with 3 rows and 9 columns
sample organism age sex infection
<character> <character> <integer> <character> <character>
GSM2545336 GSM2545336 Mus musculus 8 Female InfluenzaA
GSM2545337 GSM2545337 Mus musculus 8 Female NonInfected
GSM2545338 GSM2545338 Mus musculus 8 Female NonInfected
strain time tissue mouse
<character> <integer> <character> <integer>
GSM2545336 C57BL/6 8 Cerebellum 14
GSM2545337 C57BL/6 0 Cerebellum 9
GSM2545338 C57BL/6 0 Cerebellum 10
R
rowData(se1)
出力
DataFrame with 5 rows and 9 columns
gene ENTREZID product ensembl_gene_id
<character> <integer> <character> <character>
Asl Asl 109900 argininosuccinate ly.. ENSMUSG00000025533
Apod Apod 11815 apolipoprotein D, tr.. ENSMUSG00000022548
Cyp2d22 Cyp2d22 56448 cytochrome P450, fam.. ENSMUSG00000061740
Klk6 Klk6 19144 kallikrein related-p.. ENSMUSG00000050063
Fcrls Fcrls 80891 Fc receptor-like S, .. ENSMUSG00000015852
external_synonym chromosome_name gene_biotype phenotype_description
<character> <character> <character> <character>
Asl 2510006M18Rik 5 protein_coding abnormal circulating..
Apod NA 16 protein_coding abnormal lipid homeo..
Cyp2d22 2D22 15 protein_coding abnormal skin morpho..
Klk6 Bssp 7 protein_coding abnormal cytokine le..
Fcrls 2810439C17Rik 3 protein_coding decreased CD8-positi..
hsapiens_homolog_associated_gene_name
<character>
Asl ASL
Apod APOD
Cyp2d22 CYP2D6
Klk6 KLK6
Fcrls FCRL2
また、colData()関数を用いてサンプルメタデータから、rowData()関数を用いてフィーチャーメタデータから、それぞれ条件を指定してサブセットを抽出することもできます。例えば、ここではmiRNAと非感染サンプルのみを残します:
R
se1 <- se[rowData(se)$gene_biotype == "miRNA",
colData(se)$infection == "NonInfected"]
se1
出力
class: SummarizedExperiment
dim: 7 7
metadata(0):
assays(1): counts
rownames(7): Mir1901 Mir378a ... Mir128-1 Mir7682
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(7): GSM2545337 GSM2545338 ... GSM2545353 GSM2545354
colData names(9): sample organism ... tissue mouse
R
assay(se1)
出力
GSM2545337 GSM2545338 GSM2545343 GSM2545348 GSM2545349 GSM2545353
Mir1901 45 44 74 55 68 33
Mir378a 11 7 9 4 12 4
Mir133b 4 6 5 4 6 7
Mir30c-2 10 6 16 12 8 17
Mir149 1 2 0 0 0 0
Mir128-1 4 1 2 2 1 2
Mir7682 2 0 4 1 3 5
GSM2545354
Mir1901 60
Mir378a 8
Mir133b 3
Mir30c-2 15
Mir149 2
Mir128-1 1
Mir7682 5
R
colData(se1)
出力
DataFrame with 7 rows and 9 columns
sample organism age sex infection
<character> <character> <integer> <character> <character>
GSM2545337 GSM2545337 Mus musculus 8 Female NonInfected
GSM2545338 GSM2545338 Mus musculus 8 Female NonInfected
GSM2545343 GSM2545343 Mus musculus 8 Male NonInfected
GSM2545348 GSM2545348 Mus musculus 8 Female NonInfected
GSM2545349 GSM2545349 Mus musculus 8 Male NonInfected
GSM2545353 GSM2545353 Mus musculus 8 Female NonInfected
GSM2545354 GSM2545354 Mus musculus 8 Male NonInfected
strain time tissue mouse
<character> <integer> <character> <integer>
GSM2545337 C57BL/6 0 Cerebellum 9
GSM2545338 C57BL/6 0 Cerebellum 10
GSM2545343 C57BL/6 0 Cerebellum 11
GSM2545348 C57BL/6 0 Cerebellum 8
GSM2545349 C57BL/6 0 Cerebellum 7
GSM2545353 C57BL/6 0 Cerebellum 4
GSM2545354 C57BL/6 0 Cerebellum 2
R
rowData(se1)
出力
DataFrame with 7 rows and 9 columns
gene ENTREZID product ensembl_gene_id
<character> <integer> <character> <character>
Mir1901 Mir1901 100316686 microRNA 1901 ENSMUSG00000084565
Mir378a Mir378a 723889 microRNA 378a ENSMUSG00000105200
Mir133b Mir133b 723817 microRNA 133b ENSMUSG00000065480
Mir30c-2 Mir30c-2 723964 microRNA 30c-2 ENSMUSG00000065567
Mir149 Mir149 387167 microRNA 149 ENSMUSG00000065470
Mir128-1 Mir128-1 387147 microRNA 128-1 ENSMUSG00000065520
Mir7682 Mir7682 102466847 microRNA 7682 ENSMUSG00000106406
external_synonym chromosome_name gene_biotype phenotype_description
<character> <character> <character> <character>
Mir1901 Mirn1901 18 miRNA NA
Mir378a Mirn378 18 miRNA abnormal mitochondri..
Mir133b mir 133b 1 miRNA no abnormal phenotyp..
Mir30c-2 mir 30c-2 1 miRNA NA
Mir149 Mirn149 1 miRNA increased circulatin..
Mir128-1 Mirn128 1 miRNA no abnormal phenotyp..
Mir7682 mmu-mir-7682 1 miRNA NA
hsapiens_homolog_associated_gene_name
<character>
Mir1901 NA
Mir378a MIR378A
Mir133b MIR133B
Mir30c-2 MIR30C2
Mir149 NA
Mir128-1 MIR128-1
Mir7682 NA
課題
時刻0および時刻8におけるサンプルの先頭3つの遺伝子の発現量を抽出してください。
R
assay(se)[1:3, colData(se)$time != 4]
出力
GSM2545336 GSM2545337 GSM2545338 GSM2545341 GSM2545342 GSM2545343
Asl 1170 361 400 988 836 535
Apod 36194 10347 9173 29594 24959 13668
Cyp2d22 4060 1616 1603 3349 3122 2008
GSM2545346 GSM2545347 GSM2545348 GSM2545349 GSM2545351 GSM2545353
Asl 938 1035 494 481 937 541
Apod 27769 34301 11258 11812 29242 13682
Cyp2d22 2985 3452 1883 2014 3678 2216
GSM2545354 GSM2545380
Asl 473 1192
Apod 11088 38148
Cyp2d22 1821 4019
R
# Equivalent to
assay(se)[1:3, colData(se)$time == 0 | colData(se)$time == 8]
出力
GSM2545336 GSM2545337 GSM2545338 GSM2545341 GSM2545342 GSM2545343
Asl 1170 361 400 988 836 535
Apod 36194 10347 9173 29594 24959 13668
Cyp2d22 4060 1616 1603 3349 3122 2008
GSM2545346 GSM2545347 GSM2545348 GSM2545349 GSM2545351 GSM2545353
Asl 938 1035 494 481 937 541
Apod 27769 34301 11258 11812 29242 13682
Cyp2d22 2985 3452 1883 2014 3678 2216
GSM2545354 GSM2545380
Asl 473 1192
Apod 11088 38148
Cyp2d22 1821 4019
課題
ロング形式のrnaテーブルを使用して、同じ値が得られることを確認してください。
R
rna |>
filter(gene %in% c("Asl", "Apod", "Cyd2d22")) |>
filter(time != 4) |>
select(expression)
出力
# A tibble: 28 × 1
expression
<dbl>
1 1170
2 36194
3 361
4 10347
5 400
6 9173
7 988
8 29594
9 836
10 24959
# ℹ 18 more rows
ロングテーブルとSummarizedExperimentは同じ情報を含んでいますが、単に構造が異なるだけです。それぞれの手法に利点があります。前者はtidyverseパッケージに適しており、後者は多くのバイオインフォマティクスや統計処理の手順において好まれる構造です。例えば、DESeq2パッケージを用いた典型的なRNA-Seq解析などが挙げられます。
メタデータへの変数の追加
メタデータに情報を追加することもできます。例えば、サンプルを採取した場所を追加したいとしましょう…
R
colData(se)$center <- rep("University of Illinois", nrow(colData(se)))
colData(se)
出力
DataFrame with 22 rows and 10 columns
sample organism age sex infection
<character> <character> <integer> <character> <character>
GSM2545336 GSM2545336 Mus musculus 8 Female InfluenzaA
GSM2545337 GSM2545337 Mus musculus 8 Female NonInfected
GSM2545338 GSM2545338 Mus musculus 8 Female NonInfected
GSM2545339 GSM2545339 Mus musculus 8 Female InfluenzaA
GSM2545340 GSM2545340 Mus musculus 8 Male InfluenzaA
... ... ... ... ... ...
GSM2545353 GSM2545353 Mus musculus 8 Female NonInfected
GSM2545354 GSM2545354 Mus musculus 8 Male NonInfected
GSM2545362 GSM2545362 Mus musculus 8 Female InfluenzaA
GSM2545363 GSM2545363 Mus musculus 8 Male InfluenzaA
GSM2545380 GSM2545380 Mus musculus 8 Female InfluenzaA
strain time tissue mouse center
<character> <integer> <character> <integer> <character>
GSM2545336 C57BL/6 8 Cerebellum 14 University of Illinois
GSM2545337 C57BL/6 0 Cerebellum 9 University of Illinois
GSM2545338 C57BL/6 0 Cerebellum 10 University of Illinois
GSM2545339 C57BL/6 4 Cerebellum 15 University of Illinois
GSM2545340 C57BL/6 4 Cerebellum 18 University of Illinois
... ... ... ... ... ...
GSM2545353 C57BL/6 0 Cerebellum 4 University of Illinois
GSM2545354 C57BL/6 0 Cerebellum 2 University of Illinois
GSM2545362 C57BL/6 4 Cerebellum 20 University of Illinois
GSM2545363 C57BL/6 4 Cerebellum 12 University of Illinois
GSM2545380 C57BL/6 8 Cerebellum 19 University of Illinois
これは、メタデータスロットが他の構造に影響を与えることなく、無制限に拡張できることを示しています!
tidySummarizedExperiment
tidyverseのコマンドを使ってSummarizedExperimentオブジェクトを操作できるのか、疑問に思われるかもしれません。答えは「はい」です。tidySummarizedExperimentパッケージを使えば可能です。
SummarizedExperimentオブジェクトがどのような見た目だったか、思い出してください。
R
se
出力
class: SummarizedExperiment
dim: 1474 22
metadata(0):
assays(1): counts
rownames(1474): Asl Apod ... Lmx1a Pbx1
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(22): GSM2545336 GSM2545337 ... GSM2545363 GSM2545380
colData names(10): sample organism ... mouse center
tidySummarizedExperimentを読み込み、改めてseオブジェクトを確認してください。
R
## BiocManager::install("tidySummarizedExperiment")
library("tidySummarizedExperiment")
se
出力
class: SummarizedExperiment
dim: 1474 22
metadata(0):
assays(1): counts
rownames(1474): Asl Apod ... Lmx1a Pbx1
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(22): GSM2545336 GSM2545337 ... GSM2545363 GSM2545380
colData names(10): sample organism ... mouse center
依然としてSummarizedExperimentオブジェクトであるため、効率的な構造が維持されていますが、現在はtibbleとして扱うことができます。出力の1行目にそう記載されている点に注目してください。これはSummarizedExperiment-tibbleという抽象化です。また、出力の2行目で転写産物数とサンプル数を確認できます。
これで、tidyverseのコマンドを使ってSummarizedExperimentオブジェクトを操作できるようになりました。
filter
を使用して、条件に基づいて行をフィルタリングできます。例えば、あるサンプルのすべての行を表示する場合などです。
R
se |> filter(sample == "GSM2545336")
出力
class: SummarizedExperiment
dim: 1474 1
metadata(1): latest_filter_scope_report
assays(1): counts
rownames(1474): Asl Apod ... Lmx1a Pbx1
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(1): GSM2545336
colData names(10): sample organism ... mouse center
R
se |> filter(sex == "Female")
出力
class: SummarizedExperiment
dim: 1474 12
metadata(1): latest_filter_scope_report
assays(1): counts
rownames(1474): Asl Apod ... Lmx1a Pbx1
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(12): GSM2545336 GSM2545337 ... GSM2545362 GSM2545380
colData names(10): sample organism ... mouse center
select を使用して、表示したい列を指定できます。
R
se |> select(sample, sex)
警告
Warning: `when()` was deprecated in purrr 1.0.0.
ℹ Please use `if` instead.
ℹ The deprecated feature was likely used in the tidySummarizedExperiment
package.
Please report the issue at
<https://github.com/stemangiola/tidySummarizedExperiment/issues>.
This warning is displayed once per session.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
出力
tidySummarizedExperiment says: Key columns are missing. A data frame is returned for independent data analysis.
出力
# A tibble: 32,428 × 2
sample sex
<chr> <chr>
1 GSM2545336 Female
2 GSM2545336 Female
3 GSM2545336 Female
4 GSM2545336 Female
5 GSM2545336 Female
6 GSM2545336 Female
7 GSM2545336 Female
8 GSM2545336 Female
9 GSM2545336 Female
10 GSM2545336 Female
# ℹ 32,418 more rows
mutate を使用してメタデータ情報を追加できます。
R
se |> mutate(center = "Heidelberg University")
出力
class: SummarizedExperiment
dim: 1474 22
metadata(1): latest_mutate_scope_report
assays(1): counts
rownames(1474): Asl Apod ... Lmx1a Pbx1
rowData names(9): gene ENTREZID ... phenotype_description
hsapiens_homolog_associated_gene_name
colnames(22): GSM2545336 GSM2545337 ... GSM2545363 GSM2545380
colData names(10): sample organism ... mouse center
また、tidyverseのパイプ |>
を使ってコマンドを組み合わせることもできます。例えば、group_by
と summarise
を組み合わせることで、各サンプルの合計カウントを取得できます。
R
se |>
group_by(.sample) |>
summarise(total_counts=sum(counts))
出力
tidySummarizedExperiment says: A data frame is returned for independent data analysis.
出力
# A tibble: 22 × 2
.sample total_counts
<chr> <int>
1 GSM2545336 3039671
2 GSM2545337 2602360
3 GSM2545338 2458618
4 GSM2545339 2500082
5 GSM2545340 2479024
6 GSM2545341 2413723
7 GSM2545342 2349728
8 GSM2545343 3105652
9 GSM2545344 2524137
10 GSM2545345 2506038
# ℹ 12 more rows
整然データのSummarizedExperimentオブジェクトは、プロット用に通常のtibbleとして扱うことができます。
ここでは、サンプルごとのカウント数の分布をプロットします。
R
se |>
ggplot(aes(counts + 1, group=.sample, color=infection)) +
geom_density() +
scale_x_log10() +
theme_bw()

tidySummarizedExperimentの詳細については、パッケージのウェブサイトをご覧ください。
要点
SummarizedExperimentは、オミクスデータを保存・扱うための効率的な手法です。これらは多くのBioconductorパッケージで使用されています。
次に予定されているRNAシーケンシング解析のトレーニングでは、BioconductorのDESeq2パッケージを使用して、発現変動解析を行う方法を学びます。DESeq2パッケージを用いた解析全体は、SummarizedExperimentで行われます。
- Bioconductorは、ハイスループット生物学データの解析を支援するためのサポートやパッケージを提供するプロジェクトです。
-
SummarizedExperimentは、ハイスループットなオミクスデータを保存・管理する際に有用なオブジェクトの一種です。
Bioconductorは、R言語の2人の創始者の一人であるロバート・ジェントルマンによって立ち上げられました。Bioconductorは、オミクスデータ解析に特化したツールを提供しています。Bioconductorは統計プログラミング言語Rを採用しており、オープンソースかつオープンな開発体制となっています。↩︎