Bash on Ubuntu on WindowsでSingle-cell RNA sequence解析 ⑤番外編：コードのブラッシュアップ

※記事を読んでくださるみなさまへ

・対象はプログラミングもバイオインフォ解析も全く未経験のバイオ研究者です
・筆者は基礎研究を初めて約1年3ヶ月の大学院生です
・いわゆる”DRY”な研究を始めたのは約5ヶ月前からです
・そのため記事の内容にはとんでもない間違いやずっと良い別の方法が存在するかもしれません
・また、筆者と同じレベルの人が0から始めるときの参考になるよう、説明は結構まだるっこいです
・ご指摘の度に記事をブラッシュアップしていきたいです、よろしくお願いします

過去記事一覧
①環境構築まで
②NGSデータのダウンロードとクオリティチェック
③NGSデータのトリミング
④Pipeline (hisat2, samtools, stringtie)
⑤番外編：コードのブラッシュアップ
⑥Rによる発現解析

コードのブラッシュアップ

シェルスクリプトの勉強を進めた結果、今までの処理をもう少し自動化、簡略化できることがわかったため、今回の番外編で追記します。for文というのを利用するとよさそうです。その他前回からの変更点としては、

・出力結果のファイルを準備して、そこに書き込まれるようにした
・一つのfastqファイルをgtfファイルまで連続的に処理し、必要無くなったファイルを自動削除、ハードディスクの容量がsamファイルなどで埋め尽くされないようにした

という感じです。改行処理によりうまくいかない場合は\（環境によっては¥）を削除して一行にしてください。do, doneのところは改行をいじらないでください。

for f in *.fastq.gz
do 
java -jar ~/miniconda2/share/trimmomatic-0.36-3/trimmomatic.jar SE\
 -threads 8 ${f} ${f/.fastq.gz/_trimmed.fastq.gz}\
 ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3\
 SLIDINGWINDOW:4:15 MINLEN:36 >>result.txt\
 && fastqc ${f/.fastq.gz/_trimmed.fastq.gz} >>result.txt\
 && hisat2 -p 8 --dta -x (indexファイルのパス)/(.ht2の頭部分)\
 -U (fastqファイルのパス)/${f/.fastq.gz/_trimmed.fastq.gz}\
 -S ${f/.fastq.gz/.sam} >>result.txt\
 && samtools sort -@ 8 -o ${f/.fastq.gz/.bam}\
 ${f/.fastq.gz/.sam} >>result.txt\
 && stringtie -p 8 -G (referenceファイルのパス)/(reference.gtf)\
 -o ${f/.fastq.gz/.gtf} -l ${f/.fastq.gz/}\
 ${f/.fastq.gz/.bam} >>result.txt\
 && rm ${f}  ${f/.fastq.gz/_trimmed.fastq.gz} ${f/.fastq.gz/.sam}
done

ごちゃごちゃしていてわかりにくいですが、実は簡単です。for文は

for (変数を示すアルファベット:f,x,iなど) in (変数一覧：直接記載、'ls'、*.fastqなど)
do
(コマンド) (引数)
done

ただし、引数の中身で変数が登場するところを${f}のように記載。
変数の拡張子や名称の一部を変更したいとき${f/(変えたい部分)/(変えたいパーツ)}のようにする。

というふうな感じで、全ての候補に対して処理を行うことができます。また、&&という記号は、

(コード１) && (コード２) && ...

というふうに無限に処理をつなげて連続させることができます。また、>>は

(コード) >>〇〇.txt

でコードの出力内容をテキストファイルに追記保存していくことができます。>なら上書き保存です。

\

を入れると右隣の文字から改行をして見やすくすることができます。幾つかテストしてみたのですが、&&のような意味をもつ記号に接触する形では使わないほうが良さそうです。mkd\irみたいなのはいけます。

マージ以下は前回までと同様に処理してください。以上で番外編終了です。

（5/16追記）
どうやら今回の環境では&&をつけなくても、改行しておけば連続してコードが実行されるようです。

for f in *.fastq.gz
do 
java -jar ~/miniconda2/share/trimmomatic-0.36-3/trimmomatic.jar SE -threads 8 ${f} ${f/.fastq.gz/_trimmed.fastq.gz} ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 >>result.txt
fastqc ${f/.fastq.gz/_trimmed.fastq.gz} >>result.txt
hisat2 -p 8 --dta -x (indexファイルのパス)/(.ht2の頭部分) -U (fastqファイルのパス)/${f/.fastq.gz/_trimmed.fastq.gz} -S ${f/.fastq.gz/.sam} >>result.txt
samtools sort -@ 8 -o ${f/.fastq.gz/.bam} ${f/.fastq.gz/.sam} >>result.txt
stringtie -p 8 -G (referenceファイルのパス)/(reference.gtf) -o ${f/.fastq.gz/.gtf} -l ${f/.fastq.gz/} ${f/.fastq.gz/.bam} >>result.txt
rm ${f}  ${f/.fastq.gz/_trimmed.fastq.gz} ${f/.fastq.gz/.sam}
done