RNAseqパイプライン2.1 - kuroの覚え書き

フルパス付きとか拡張子とか気にせずにすべてドラッグ＆ドロップでできるようにスクリプトをブラッシュアップ
tophat_cufflinks2

#!/bin/sh

#アノテーション情報を含むgtfファイル
read -p "annotation infomation (full path) = " gtf
#リファレンスファイルの指定
read -p "reference file (full path) = " ref
echo "アノテーション情報は"$gtf""
echo "リファレンスファイルは"$ref""
for file in $@
do
#ファイル名を取り出す
fname="${file##*/}"
echo ""$file"についてtophatでmappingを実行する"
#出力するディレクトリを作る
dir_name="${fname%.*}"
mkdir -p ./tophat_results/"$dir_name"
output_dir=./tophat_results/"$dir_name"
#echo "出力先は./tophat_results/"$dir_name""
tophat -p 4 -G "$gtf" -o "$output_dir" "${ref%.*}" "$file"
echo "---------------------------------------------"
echo "accepted_hits.bamを"$output_dir"/"$dir_name"_th.bamに改名する"
mv "$output_dir"/accepted_hits.bam "$output_dir"/"$dir_name"_th.bam
echo "---------------------------------------------"
echo ""$dir_name"_th.bamについてcufflinksを実行する"
output_dir2=./tophat_results/"$dir_name"/cufflinks_results
cufflinks -p 4 -g "$gtf" "$output_dir"/"$dir_name"_th.bam -o "$output_dir2"
echo "---------------------------------------------"
done

使用方法

$ cd /Volumes/HDD1/expression/hogehoge/trimmed/
#トリミングできたfastqファイルをおいているディレクトリに移動(任意のディレクトリで構わないがその下に結果ファイルができる）
$ tophat_cufflinks2 /Volumes/HDD1/expression/hogehoge/trimmed/hoge1_trim.fastq /Volumes/HDD1/expression/hogehoge/trimmed/hoge2_trim.fastq /Volumes/HDD1/expression/hogehoge/trimmed/hoge3_trim.fastq

ファイルはターミナルウィンドウにドラッグ＆ドロップで追加していってリターンすると・・・

annotation infomation (full path) = /Users/hoge/expression/Homo_sapiens/NCBI/build37.2/Annotation/Archives/archive-2014-06-02-13-47-29/Genes/genes.gtf 
reference file (full path) = /Users/hoge/expression/Homo_sapiens/NCBI/build37.2/Sequence/Bowtie2Index/genome.fa

アノテーション情報ファイルとリファレンスファイルをこれまたドラッグ＆ドロップで指定

アノテーション情報は/Users/hoge/expression/Homo_sapiens/NCBI/build37.2/Annotation/Archives/archive-2014-06-02-13-47-29/Genes/genes.gtf
リファレンスファイルは/Users/hoge/expression/Homo_sapiens/NCBI/build37.2/Sequence/Bowtie2Index/genome.fa
/Volumes/WD1/expression/ERR266a/trimmed/ERR266335a_trim.fastqについてtophatでmappingを実行する

[2017-05-10 10:41:11] Beginning TopHat run (v2.1.1)
-----------------------------------------------
[2017-05-10 10:41:11] Checking for Bowtie
		  Bowtie version:	 2.3.0.0

以下略

こんな感じに進んでいく

同様にトリミングもD&Dで
fastqc_trim_qc2

#$!/bin/sh
#
#QC用データ格納ディレクトリがなければ作成し、QCを実施
#クオリティでフィルター後、トリミング（パラメータは本スクリプトを直接書き換えておく）
#QCを再度かける
#GNU paralleleを用いて並列処理するには
# $ parallele fastqc_trim_qc2 ::: hoge1 hoge2 hoge3 (ファイル名、フルパスも可）
#
mkdir -p ./FastQC/trimmed
mkdir -p ./trimmed
#Phred quality score "qscore"未満がconc%以上のリードを除去
#Phred quality score "qscore"未満の末端を除去
    qscore=20
    conc=80
#配列長が"length"bp未満のリードを除去
    length=30
echo "minimum Phred quality score = "$qscore""
echo "minimum base length = "$length""
echo "% threshfold = "$conc""
for file in "$@"
do
#ファイル名を取り出す
fname="${file##*/}"
name_pref="${fname%.*}"
#FastQCにかける
fastqc --nogroup -o ./FastQC "$fname"
fastq_quality_filter -Q33 -q "$qscore" -p "$conc"  -i "$fname" | fastq_quality_trimmer -Q33 -t "$qscore" -l "$length" -o trimmed/"$name_pref"_trim.fastq
#FastQCにかける
fastqc --nogroup -o ./FastQC/trimmed trimmed/"$name_pref"_trim.fastq
done