前処理

パッケージのインストール

conda install -c bioconda trimmomatic
conda install -c bioconda fastqc

アダプター配列のトリミングとQC

trimmomaticを用いてトリミングを行う。
trimmomaticでは、アダプター配列を指定でき、いくつかの配列はプリセットで保存されている。
アダプター配列ファイルへのシンボリックリンクを作成し、カレントディレクトリ上のファイルとして参照可能とする。
ln -s /Users/minami/opt/anaconda3/share/trimmomatic-0.39-2/adapters/TruSeq3-SE.fa

trimmomatic.

FILES=(今回使用した4種の条件名を記載。半角スペースをあけて表記する)
for file in ${FILES[@]}; do
trimmomatic PE -threads 10 -phred33 -trimlog trimlog_${file}.txt -summary summary_${file}.txt ${file}_1.fq.gz ${file}_2.fq.gz ${file}_trim_pair_1.fastq.gz ${file}_trim_unpair_1.fastq.gz ${file}_trim_pair_2.fastq.gz ${file}_trim_unpair_2.fastq.gz ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:30
done

fastqc.

mkdir fastqc
for file in ${FILES[@]}; do
fastqc -t 10 --nogroup -o fastqc -f fastq {file}_trim_pair_1.fastq.gz ${file}_trim_pair_2.fastq.gz
done

マッピングとソーティング

パッケージのインストール

conda install -c bioconda hisat2
conda install -c bioconda samtools

samtoolsのインストール時の問題

samtools
とすると、以下のエラーが出た。
samtools: error while loading shared libraries: libcrypto.so.1.0.0: cannot open shared object file: No such file or directory
こちらの記事を参考に、シンボリックリンクを作成した。

~/seq$ cd /home/member/anaconda3/bin/../lib
~/anaconda3/lib$ ln -s libcrypto.so.1.1 libcrypto.so.1.0.0

hisat2を用いたマッピング

手持ちのMacbook Air(M1, メモリ8GB)ではまったく進まなかった。メモリ不足か。ラボ共用の高スペックPCで行ったところすぐに実行できた。

for file in ${FILES[@]};do
> hisat2 -p 16 -t -x grcm38/genome -1 ${file}_trim_pair_1.fastq.gz -2 ${file}_trim_pair_2.fastq.gz -S ${file}.sam
> done

マッピング結果の確認（一例）

22381736 reads; of these:
  22381736 (100.00%) were paired; of these:
    1221246 (5.46%) aligned concordantly 0 times
    17904164 (79.99%) aligned concordantly exactly 1 time
    3256326 (14.55%) aligned concordantly >1 times
    ----
    1221246 pairs aligned concordantly 0 times; of these:
      135217 (11.07%) aligned discordantly 1 time
    ----
    1086029 pairs aligned 0 times concordantly or discordantly; of these:
      2172058 mates make up the pairs; of these:
        1254782 (57.77%) aligned 0 times
        740338 (34.08%) aligned exactly 1 time
        176938 (8.15%) aligned >1 times
97.20% overall alignment rate  

となった。90％超えていれば問題なさそう？

samtoolsを用いてsamファイルをソート済bamファイルに変換

for file in ${FILES[@]};do 
samtools sort -@ 16 -O bam -o ${file}_sort.bam ${file}.sam; 
done

ソート済bamファイルのindexを作成

indexが必要になるツールもあるらしい（まだやってない）

for file in ${FILES[@]};do samtools index ${file}_sort.bam; done

stringtieを用いた発現量の推定

あとでballgownを用いて解析するので、ballgown用のオプションを指定しておく

for file in ${FILES[@]};do 
stringtie ${file}_sort.bam -e -B -G gencode.vM28.annotation.gtf -o ballgown/${file}/${file}.gtf
done

ballgownを用いた発現変動遺伝子の抽出

ここからはRを用いる

#biocmanagerからballgownパッケージをインストール
install.packages("BiocManager")
BiocManager::install("ballgown")

#途中でUpdate all/some/none?と聞かれたら、「a」を入力してenter

install.packages('metaMA')
library(metaMA)

install.packages('tidyr')
library(tidyr)

install.packages('dplyr')
library(dplyr)

library(ballgown)

#ballgownへのデータの読み込み
#公的データベースからダウンロードしたデータの場合基本的に頭文字3文字が同じだが、委託解析したデータでは違う場合がある
#そのときはballgownフォルダ内の名前に頭文字3文字を適当につければ良い
#samplePatternにはその文字を指定する

b <- ballgown(dataDir = 'ballgown', samplePattern = 'SRR')

#id,groupは適当なものを指定
pData(b) <- data.frame(id =  c('ctrl_1', 'ctrl_2', 'ctrl_3', 'mut_1', 'mut_2', 'mut_3'),
                       group = c('ctrl', 'ctrl', 'ctrl', 'mut', 'mut', 'mut'))

転写産物の発現量FPKMを取得したい場合はtexpr関数を、遺伝子ごとの発現量を取得する場合はgexpr関数を利用する。 発現量の分散がある程度大きい遺伝子を抽出する。

b_highexp <- subset(b, 'rowVars(texpr(b)) > 1', genomesubset = TRUE)
tx_highexp <- texpr(b_highexp)

ここからstattest関数を用いると、発現変動遺伝子を抽出できる。
今回はn=1で行ったため、stattestではなく、tx_highexpをcsvで出力し、pythonで処理することで変動遺伝子候補を抽出した。