Python学んで人生変わった（バイオ系博士)
今回の内容
私は、学部生の頃はスポーツ科学を学んでおり、大学院から専攻を基礎生物学に変えました。パソコン苦手の文系の私が、pythonを学び、就職活動を終えるまでをまとめてみました。
大学院に入るまでの私の能力
- 大学時代まで体育学部（脳みそ筋肉じゃないと信じたい笑）
- プログラミング言語なんて聞いたこともなかった
- パソコン苦手のバリバリの文系
始まりの日
「じゃあ、プログラミングを使って解析していこうか」
プログラミング書けるって良さげじゃねという軽い気持ちで始めた最初の日。
本当に、プログラミング言語って何というところから始まりました（笑）
読めない、書けない、わからない。
ですが、研究を進めるためには、生物系の解析でよく使われるPython（プログラミング言語）を習得しなければなりません。当時の私はモチベーションだけはありました。。。
待ち受ける困難
さっぱりわからない。本当にわからない。print構文しかわからない (笑)
Pythonは比較的、簡単と言われるのに全然できない（汗）
こんなの一生できなよと思った日々orz
最高潮にあったモチベーションは徐々に落ちはじめました。。。
ですが、プログラミングが書けなければ研究ができません。なので、まずscriptを読む練習から始めました (まるで、英文読解)。
思わぬ幸運
エラーばかり出て発狂しそうになったこともありましたが、ネット上にあるスクリプトを使って読む練習をしていくと、なぜか少しずつわかるようになってきました（笑）
順序としては、
- あ、この一文わかるぞ
- だいぶ読めてきたぞ、しかし、書けない。。。
- これなら書けそう！ (ほんとたまに)
- おぉー書ける、書けるぞー！
そして、Pythonを学んで、知らず知らずのうちに得たこととは、
そう、就活の幅が広がる
そして、バイオ系の企業にも入りやすくなる
これは、とても大きかったです。特にバイオ系の研究していると就職が難しいと言う話をよく聞きますが、プログラミングができるだけで、かなり違います。
最近ですと、バイオ系の企業もプログラミングを書ける人材を募集していますし、システムエンジニア(SE)など情報系の職にも応募できるようになります
新たなる旅立ち
修士、博士の５年間で生物の研究を行い、もう社会人では生物系の研究はせず、一般職を選ぼうと思っていました。
そんな私が、企業の研究職に着き、今後も生物系のバイオインフォマティクス解析を行うことになりました（笑）
シングルセル解析におけるCSVファイルの取り込み（scanpy）
シングルセル解析を行う際、scanpyを使った解析をよく見ます。その中で、CSVファイルを取り込む方法が、情報として少なかったため、まとめて見ました。
通常の解析に関しては、【Python】Scanpyを使った single cell RNA解析 - ばいばいバイオに詳しく書かれています。
Scanpyとは？
Scanpy is a scalable toolkit for analyzing single-cell gene expression data built jointly with anndata.
シングルセル用のデータ解析ツールであり、anndata形式に対応している。
Scanpyより引用
Adataとは？
デモデータを参考にすると、
var_namesには遺伝子名が入ります。
.obsには細胞名やクラスタリング後のクラスター番号が加えられます。
.Xには、数値データのみが入り、マトリックス形式となっています。
scanpyを用いて解析途中のデータを保存する際は、拡張子をh5adにし、保存する方がその後の解析がしやすそうです。
COVID-19のdatasetsもh5adらしい 2020/09/26現在
/ adata.write(file_path, compression='gzip') #ファイルの拡張子を変えて、ファイル形式を指定できる compressionは圧縮形式を指定 /
CSVファイルの取り込み
adata = sc.read(file_path, cache =True) #loading csv faster than pd.csv adata.obs=adata.obs.reset_index(drop=True) #resetting adata.obs なぜかエラーが出てくる cell_type=pd.read_csv(File_Path[2],delimiter="\t") adata.obs['cell_type']=cell typeのリスト #obsにcell typeを入れる gene_id=pd.read_csv(File_Path[3],delimiter="\t") adata.var['gene_id']=Gene_idのリスト #varにgene idを入れる
参考
StringTieを用いたリード数の計算 (RNAseq)
マッピング後のSAM、BAMから、各遺伝子ごとにリード数を算出します。今回は、StringTieについてまとめてみました。
SAMファイルは重たいので、BMAファイルに変換した方が解析がしやすいと思います。
The main input of the program is a BAM file with RNA-Seq read mappings which must be sorted by their genomic location (for example the accepted_hits.bam file produced by TopHat or the output of HISAT2
準備
マッピング後のファイルは事前にソートを行なっておく
samtools sort - input.sorted
実行コマンド
stringtie [input_sort.bam] [options]
使いそうな設定をいくつか抜粋
| オプション | 説明 |
|---|---|
| -o | 出力するGTFファイルパスの設定 |
| -p | スレッド数 デフォルト１ |
| -G | GTF、GFF3 |
| -f | 最も多く存在するアイソフォームに対して、最小のアイソフォーム転写物の存在量を設定 デフォルト0.01 |
| -m | 新規転写物の最小長を設定 デフォルトは２００ |
| -c | 出力の転写物名をレファレンスを基に表示する -Gオプションが必要 |
| -s | エキソンが一つしかない転写物の最小カバレッジ量を指定する |
| -A | タブ区切りの出力を行う際の保存場所の指定（出力形式は下記を参照） |
GTF、GFFファイルを指定する場合は、マッピングに使用したゲノムの染色体名と一致しているかを確認する。
小文字、大文字の違いで上手く処理できない場合がある。
例文
stringtie -p 8 -G TAIR10_for_stringtie.gff -o count.txt -A table.txt sorted.bam
出力形式
- seqname: マップされたクロモソーム名
- source: 基となったGTFファイル
- feature: RNAの特徴 (e.g., exon, transcript, mRNA, 5'UTR).
- start: 開始位置 1-based index
- end: 終了位置 1-based index
- score: 現在は使用されていない。デフォルト表記は１０００？
- strand: ストランド
- frame: 表記なし（本来はCDSのフレーム）
- attributes: その他の特徴を記載
FPKM、TPMに関しては下記のサイトを参照
RNA-seqブログ
-Aを用いた出力
- Column 1 / 遺伝子ID
- Column 2 / 遺伝子名
- Column 3 / 使用したReferenceの名前
- Column 4 / ストランド
- Column 5 / 開始位置 (1-based index).
- Column 6 / 終了位置 (1-based index).
- Column 7 / Coverage 各塩基単位のリード数
- Column 8 / FPKM
- Column 9 / TPM
参考
SRAからシーケンスデータをダウンロードする（RNAseq）
NCBI Sequence Read Archive (SRA) とは？
Sequence Read Archive (SRA) data, available through multiple cloud providers and NCBI servers, is the largest publicly available repository of high throughput sequencing data. SRA stores raw sequencing data and alignment information to enhance reproducibility and facilitate new discoveries through data analysis.
Sequence Read Archive (SRA) データは、複数のクラウドプロバイダー、NCBI サーバーを介して利用可能で、一般に公開されている中でハイスループットシーケンシングデータを扱う最大のリポジトリです。SRAは、生シーケンシングデータとアラインメント情報を保存することで、再現性を高め、データ解析による新たな発見を促進します。
NCBIより引用
sratoolkitのインストール
brew install sratoolkit
ファイルのダウンロード、変換
prefetch SRR number (e.g. SRR7652713) シングルエンド（SE） fastq-dump --gzip directory_of_SRR ペアエンド（PE） fastq-dump --gzip --split-files directory_of_SRR
参考
SAMファイルの変換（RNAseq）
マッピング後の処理（SAM、BAM）
samファイルはシーケンスリードがゲノムのどの位置にあるかを記述するファイル形式であり、様々なシーケンス解析のツールで利用されている。samファイルはテキスト形式のデータであり、ファイルサイズが大きくなったり、ソフトウェアから効率的に個々のデータにアクセスすることが難しい。そこでbam形式というバイナリ形式のファイルを作成することで、ファイルサイズの縮小、アクセスの高速化、リードのゲノム座標でのソートによるアクセスの効率化が期待できる。bai形式は bam のインデックスファイルで、このファイルをプログラムが参照することで、効率的に bam file にアクセスできる。
Tips for NGS Data Analysis より引用
マッピング後、SAMファイルができます。このSAMファイルをBAMファイルに変換することで、ファイルサイズを縮小し、cufflinks、stringtieなどを用いてRNA蓄積量を算出することができます。
また、BAMファイルをBEDファイルに変換することでマッピングされたリードの末端情報を取得することができます。
加えて、BAMファイルをソートし、インデックス処理を行うことで、マッピングされたリードをIGVなどを用いて視覚化することができます。
SAMファイルをBAMファイルに変換
samtools view -Sb file.sam > file.bam
BAMファイルをBEDファイルに変換
bedtools bamtobed -i file.bam > file.bed
BAMファイルをソート、インデックスを付加
samtools sort -@ 8 -o file.sort file.bam #bamファイルをソート or samtools sort -@ 8 file.bam -o file #bamファイルをソート(上書きする) samtools index file.sort.bam #インデックスを作成 bam.baiを作成
参考
Hisat2を用いたマッピング (RNAseq)
indexの作成
hisat2-build [解析対象とするゲノムfastaファイル] [保存先を指定] e.g. hisat2-build Arabidopsis_thaliana.TAIR10.dna.fa TAIR10
基本設定
シングルエンド（SE） hisat2 [options] -x [作成したindexファイル] -U [解析対象とするfastaファイル] -S [出力先を指定] ペアエンド（PE） hisat2 [options] -x [作成したindexファイル] -1 [リード１] -2 [リード２] -S [出力先を指定]
詳細な設定
option（共通） 説明 -p スレッド数の変更 デフォルト１ -s スキップするリード数を指定 -u 指定した数をinputとして使用 -5 5' 末端から指定した数を削る -3 3' 末端から指定した数を削る --no-softclip ソフトクリップを解析から除外 --score-min 最小のアライメントスコアを指定する --min-intronlen イントロンの最小長を指定 --max-intronlen イントロンの最長を指定 --known-splicesite-infile 既知のスプライシングサイトを考慮する --novel-splicesite-outfile スプライシングサイトを出力する --no-unal マッピングされなかったリードを出力しない --un-gz アンマップリードを出力する --summary-file サマリーを出す --new-summary 新形式のサマリーを出す
PE 説明 -I ペアエンドでマップされた際の最小長を指定 デフォルト0 -X ペアエンドでマップされた際の最大長を指定 デフォルト５００ --no-mixed ペアでマップされなかった場合、各々でアライメント箇所の探索を行わない --no-discordant リード１とリード２が適切な順序、領域でない場合（遺伝子をまたぐなど）、アライメントを行わない
参考
Fastx tool kitを用いたフィルタリング (RNAseq)
個人的にはcutadaptやtrimomaticが使いやすいと思いますが、先行研究などと合わせたい場合に使用することがあるため、Fastx tool kitに関してまとめてみました！
[-h] 使い方を表示
[INFILE] input ファイルを指定 (FASTAファイルの場合、塩基比率のみ算出)
[OUTFILE] アウトプットをテキストファイルで出力
fastx_quality_stats [-h] -i [INFILE] -o [OUTFILE]
fastx_nucleotide_distribution_graph.sh -i [input file] -t [title name] [-p] -o [output file]
fastx_clipper [-h] -a [ADAPTER] -I [min_length] -d [N] [-k] [-c] [-C] [-n] [-v] [-z] -i [INFILE] -o [OUTFILE]
[-h] ヘルプを表示 [ADAPTER] アダプター配列を指定 デフォルトはCCTTAAGG
[min_length] 指定した塩基長より短い場合、解析から除外 デフォルトは5
[N] アダプター配列以降の配列をN数塩基保持？
[-c] アダプター配列が存在しないリードを解析から除外
[-C] アダプター配列が存在したリードを解析から除外
[-k] アダプター配列だけのリードを報告する？
[-n] Nが存在する配列を保持 デフォルトはNが持つリードを除去
[-v] リード数を報告
[-z] outputファイルをGZIP圧縮する
[INFILE] input ファイルを指定
[OUTFILE] output ファイルを指定
fastx_trimmer [-h] -f [first_base] -I [last_base] [-z] [-v] -i [INFILE] -o [OUTFILE]
[-h] ヘルプを表示
[first_base] 指定した塩基以降を保持する デフォルトは1 (１塩基目)
[last_base] 指定したlast_base以降の塩基を除去 デフォルトは全塩基を保持
[-z] GZIP で圧縮
[INFILE] input ファイルを指定
[OUTFILE] output ファイルを指定
fastq_quality_filter [-h] [-v] -q [min_quality] -p [min_quality_percent] [-z] -i [INFILE] -o [OUTFILE]
[-h] ヘルプを表示
[min_quality] N以下のクオリティ塩基を除去する
[min_quality_percent] [min_quality] 以下のクオリティ塩基がN％以上のリードを除去
[-z] GZIP形式で圧縮
[INFILE] input ファイルを指定
[OUTFILE] output ファイルを指定
[-v] リード数を報告する
