Merge pull request #33 from yyoshiaki/developmentv1.2

v1.2
yyoshiaki · May 26, 2019 · 38e2449 · 38e2449
2 parents b0f6589 + 39b2307
commit 38e2449
Show file tree

Hide file tree

Showing 4 changed files with 100 additions and 28 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/README.md b/README.md
@@ -1,4 +1,5 @@
-# ikra v1.1.1 -RNAseq pipeline centered on Salmon-<img src="img/ikra.png" width="20%" align="right" />
+# ikra v1.2.0 -RNAseq pipeline centered on Salmon-<img src="img/ikra.png" width="20%" align="right" />
+
 
 [idep](http://bioinformatics.sdstate.edu/idep/)のinputとして発現量テーブル（gene × sample）をexperiment matrixから自動でつくる。salmonを用いる。
 
@@ -9,8 +10,8 @@ ikraの`tximport_R.R`にサンプルを取り違えうる重大なバグが見
 ## Usage
 
 ```
-Usage: ikra.sh experiment_table.csv spiece \
-        [--test, --fastq, --help, --without-docker, --udocker] \
+Usage: ikra.sh experiment_table.csv species \
+        [--test, --fastq, --help, --without-docker, --udocker --protein-coding] \
         [--threads [VALUE]][--output [VALUE]]\
         [--suffix_PE_1 [VALUE]][--suffix_PE_2 [VALUE]]
   args
@@ -22,11 +23,13 @@ Options:
   --fastq use fastq files instead of SRRid. The extension must be foo.fastq.gz (default : False)
   -u, --udocker
   -w, --without-docker
+  -pc, --protein-coding use protein coding transcripts instead of comprehensive transcripts.
   -t, --threads
   -o, --output  output file. (default : output.tsv)
   -s1, --suffix_PE_1    suffix for PE fastq files. (default : _1.fastq.gz)
   -s2, --suffix_PE_2    suffix for PE fastq files. (default : _2.fastq.gz)
   -h, --help    Show usage.
+  -v, --version Show version.
 ```
 
 1. test optionは各サンプルにおいてリード数を100000に限定する。
@@ -53,27 +56,33 @@ experiment matrixはカンマ区切りで（csv形式）。
 
 - nameはアンダーバー区切りでcondition、replicateをつなげて書く。
 - 前3列は必須。
-- 自前のfastq fileを使いたいときは`--fastq`をつける。拡張子は`fastq.gz`のみに対応。
+- 自前のfastq fileを使いたいときは`--fastq`をつける。拡張子は`.fq`, `.fq.gz`, `.fastq`, `fastq.gz`のみに対応。
 - fastq fileは`fastq.gz`もしくは`_1.fastq.gz`,`_2.fastq.gz`を除いたpathを。例えば`hoge/SRR5385247.fastq.gz`なら`hoge/SRR5385247`と記載。
 - suffixが`_1.fastq.gz`,`_2.fastq.gz`ではない場合は-s1, -s2オプションをつける。
+- `../fq/**.fastq.gz`など、実行ディレクトリより上の階層を指定することはdockerの都合上不可能だが、symbolic linkを貼ることで回避できる。
+[bonohu blog](https://bonohu.github.io/running-ikra.html)
 
 - Illumina用 : trimmomatic -> trim_galoreに切り替えた。
 - Ion S5用: SEしか無い。trimmomaticではなくfastx-toolsを使う。adapterはNoneを入れておく。(test : [DRP003376](https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=DRP003376))
 
 ### Output
 
-- output.tsv
+- output.tsv(scaledTPM)
 
 - multiqc_report.html
 salmonのマッピング率（トランスクリプトに対するマッピング率）
 
-### 仕様について
+### 各種仕様
 
 - outputは**scaledTPM** (see. [Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4, 1521 (2015).](https://f1000research.com/articles/4-1521/v2))。
 - GCbiasについて、salmonで`--gcBias`を追加した。GCbiasのRNAseqにおける影響に関しては[Mike Love's blog :
 RNA-seq fragment sequence bias](https://mikelove.wordpress.com/2016/09/26/rna-seq-fragment-sequence-bias/)。
 - validateMappings optionを採用。（alignment-base modeでは使えない。）詳しくは[salmon Frequently Asked Questions](https://combine-lab.github.io/salmon/faq/)。
 
+## 重要　bugについて　2019/04/30
+
+ikraの`tximport_R.R`にサンプルを取り違えうる重大なバグが見つかり、修正しました。必ずv1.1.1以降に更新してお使いください。古いバージョンを使われていた方は、中間ファイルは問題ありませんので、`output.tsv`を削除し、もう一度新しいikra.shを実行してください。大変ご迷惑をおかけいたしました。
+
 ## Install
 
 dockerかudocker(v1.1.3)をインストール済みであること。
@@ -109,7 +118,7 @@ $ cd test/Illumina_SE && bash ../../ikra.sh Illumina_SE_fastq.csv mouse --fastq
 **SRR mode**
 
 ```bash
-$ cd test/Illumina_PE && bash ../../ikra.sh Illumina_PE_SRR.csv mouse --test -t 50
+$ cd test/Illumina_PE && bash ../../ikra.sh Illumina_PE_SRR.csv mouse --test -t 10
 ```
 
 **fastq mode**
@@ -149,6 +158,8 @@ SRRデータを探している場合は[http://sra.dbcls.jp/](http://sra.dbcls.j
 
 ## やったこと
 
+詳しくは[Relases](https://github.com/yyoshiaki/ikra/releases)を参照。
+
 - udockerの対応
 - 生物種の判別(アナログ)
 - gtf, transcript file をGENCODEから
@@ -164,6 +175,7 @@ SRRデータを探している場合は[http://sra.dbcls.jp/](http://sra.dbcls.j
 - fasterq-dump
 - cwl開発少しだけ
 - 名前の変更（ikra）
+- protein coding option
 
 ## legacy
 
@@ -191,7 +203,6 @@ trimmomaticを使ったトリミングを用いたフローは`./legacy`に移
 cd test/cwl_PE && bash test.sh
 ```
 
-
 ## cwl_toolsの由来、参考
 
 - https://github.com/pitagora-galaxy/cwl

diff --git a/ikra.sh b/ikra.sh
@@ -9,7 +9,7 @@ COMMENTOUT
 #　オプション関連ここから
 #　大部分は http://dojineko.hateblo.jp/entry/2016/06/30/225113 から引用させていただきました。
 
-#　変数 EX_MATRIX_FILE, REF_SPIECE はここで定義
+#　変数 EX_MATRIX_FILE, REF_SPECIES はここで定義
 #　if [[ $IF_TEST = true ]]; then でテストモード用の実行が可能
 
 #　今まで$1 = EX_MATRIX_FILEだったのを変更している
@@ -18,11 +18,14 @@ COMMENTOUT
 set +u
 
 PROGNAME="$( basename $0 )"
+VERSION="v1.2.0dev"
 
 # Usage
 function usage() {
   cat << EOS >&2
-Usage: ${PROGNAME} experiment_table.csv spiece [--test, --fastq, --help, --without-docker, --udocker] [--threads [VALUE]][--output [VALUE]][--suffix_PE_1 [VALUE]][--suffix_PE_2 [VALUE]]
+ikra ${VERSION} -RNAseq pipeline centered on Salmon
+
+Usage: ${PROGNAME} experiment_table.csv species [--test, --fastq, --help, --without-docker, --udocker, --protein-coding] [--threads [VALUE]][--output [VALUE]][--suffix_PE_1 [VALUE]][--suffix_PE_2 [VALUE]]
   args
     1.experiment matrix(csv)
     2.reference(human or mouse)
@@ -32,11 +35,22 @@ Options:
   --fastq use fastq files instead of SRRid. The extension must be foo.fastq.gz (default : False)
   -u, --udocker
   -w, --without-docker
+  -pc, --protein-coding use protein coding transcripts instead of comprehensive transcripts.
   -t, --threads
   -o, --output  output file. (default : output.tsv)
   -s1, --suffix_PE_1    suffix for PE fastq files. (default : _1.fastq.gz)
   -s2, --suffix_PE_2    suffix for PE fastq files. (default : _2.fastq.gz)
   -h, --help    Show usage.
+  -v, --version Show version.
+EOS
+  exit 1
+}
+
+
+# version
+function version() {
+  cat << EOS >&2
+ikra ${VERSION} -RNAseq pipeline centered on Salmon
 EOS
   exit 1
 }
@@ -47,6 +61,7 @@ DOCKER=docker
 THREADS=1
 IF_TEST=false
 IF_FASTQ=false
+IF_PC=false
 SUFFIX_PE_1=_1.fastq.gz
 SUFFIX_PE_2=_2.fastq.gz
 OUTPUT_FILE=output.tsv
@@ -62,6 +77,9 @@ for opt in "$@"; do
         '--fastq' )
             IF_FASTQ=true; shift
             ;;
+        '-pc'|'--protein-coding' )
+            IF_PC=true; shift
+            ;;
         '-u'|'--undocker' )
             DOCKER=udocker; shift
             ;;
@@ -104,6 +122,9 @@ for opt in "$@"; do
         '-h' | '--help' )
             usage
             ;;
+        '-v' | '--version' )
+            usage
+            ;;
         '--' | '-' )
             shift
             PARAM+=( "$@" )
@@ -123,10 +144,10 @@ done
 
 # オプション無しの値を使う場合はここで処理する
 EX_MATRIX_FILE="${PARAM}"; PARAM=("${PARAM[@]:1}")
-REF_SPIECE="${PARAM}"; PARAM=("${PARAM[@]:1}")
+REF_SPECIES="${PARAM}"; PARAM=("${PARAM[@]:1}")
 
 [[ -z "${EX_MATRIX_FILE}" ]] && usage
-[[ -z "${REF_SPIECE}" ]] && usage
+[[ -z "${REF_SPECIES}" ]] && usage
 
 # 規定外のオプションがある場合にはusageを表示
 if [[ -n "${PARAM[@]}" ]]; then
@@ -136,12 +157,13 @@ fi
 # 結果を表示(オプションテスト用)
 cat << EOS | column -t
 EX_MATRIX_FILE ${EX_MATRIX_FILE}
-REF_SPIECE ${REF_SPIECE}
+REF_SPECIES ${REF_SPECIES}
 RUNINDOCKER ${RUNINDOCKER}
 DOCKER ${DOCKER}
 THREADS ${THREADS}
 IF_TEST ${IF_TEST:-false}
 IF_FASTQ ${IF_FASTQ:-false}
+IF_PC ${IF_PC:-false}
 EOS
 
 set -u
@@ -156,21 +178,31 @@ SRA_ROOT=$HOME/ncbi/public/sra
 
 SCRIPT_DIR=$(cd $(dirname $0); pwd)
 
-if [[ $REF_SPIECE = mouse ]]; then
-  BASE_REF_TRANSCRIPT=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M19
-  REF_TRANSCRIPT=gencode.vM19.transcripts.fa.gz
+if [[ $REF_SPECIES = mouse ]]; then
+  BASE_REF_TRANSCRIPT=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M21
+  REF_TRANSCRIPT=gencode.vM21.transcripts.fa.gz
+  if [ $IF_PC = false ]; then
+    REF_TRANSCRIPT=gencode.vM21.transcripts.fa.gz
+  else
+    REF_TRANSCRIPT=gencode.vM21.pc_transcripts.fa.gz
+  fi
   SALMON_INDEX=salmon_index_mouse
-#   REF_GTF=gencode.vM19.annotation.gtf.gz
-  TX2SYMBOL=gencode.vM19.metadata.MGI.gz
+#   REF_GTF=gencode.vM21.annotation.gtf.gz
+  TX2SYMBOL=gencode.vM21.metadata.MGI.gz
+
+elif [[ $REF_SPECIES = human ]]; then
+  BASE_REF_TRANSCRIPT=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30
+  # REF_TRANSCRIPT=gencode.v30.pc_transcripts.fa.gz
+
+  if [ $IF_PC = false ]; then
+    REF_TRANSCRIPT=gencode.v30.transcripts.fa.gz
+  else
+    REF_TRANSCRIPT=gencode.v30.pc_transcripts.fa.gz
+  fi
 
-elif [[ $REF_SPIECE = human ]]; then
-  BASE_REF_TRANSCRIPT=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29
-  # REF_TRANSCRIPT=gencode.v29.pc_translations.fa.gz
-  REF_TRANSCRIPT=gencode.v29.transcripts.fa.gz
   SALMON_INDEX=salmon_index_human
 #   REF_GTF=gencode.v29.annotation.gtf.gz
-  TX2SYMBOL=gencode.v29.metadata.HGNC.gz
-
+  TX2SYMBOL=gencode.v30.metadata.HGNC.gz
 else
   echo No reference speice!
   exit
@@ -404,6 +436,16 @@ fi
 # trim_galore
 # SE
 if [ $LAYOUT = SE ]; then
+  if [[  -f "${dirname_fq}${SRR}.fq" ]]; then
+    mv ${dirname_fq}${SRR}.fq ${dirname_fq}${SRR}.fastq
+  fi
+  if [[  -f "${dirname_fq}${SRR}.fastq" ]]; then
+    $PIGZ ${dirname_fq}${SRR}.fastq.gz
+  fi
+  if [[  -f "${dirname_fq}${SRR}.fq.gz" ]]; then
+    mv ${dirname_fq}${SRR}.fq.gz ${dirname_fq}${SRR}.fastq.gz
+  fi
+
   if [[ ! -f "${dirname_fq}${SRR}_trimmed.fq.gz" ]]; then
     $TRIMGALORE ${dirname_fq}${SRR}.fastq.gz
   fi
@@ -416,7 +458,7 @@ if [ $LAYOUT = SE ]; then
 # PE
 else
   # trimmomatic
-  if [[ ! -f " ${dirname_fq}${SRR}_1_val_1.fq.gz" ]]; then
+  if [[ ! -f "${dirname_fq}${SRR}_1_val_1.fq.gz" ]]; then
     $TRIMGALORE --paired ${dirname_fq}${SRR}${SUFFIX_PE_1} ${dirname_fq}${SRR}${SUFFIX_PE_2}
   fi
 
@@ -513,6 +555,10 @@ if [[ ! -f "$OUTPUT_FILE" ]]; then
   $RSCRIPT_TXIMPORT tximport_R.R $TX2SYMBOL $EX_MATRIX_FILE $OUTPUT_FILE
 fi
 
+# tximport
+if [[  -f "tximport_R.R" ]]; then
+  rm tximport_R.R
+fi
 
 # if [[ "$RUNINDOCKER" -eq "1" ]]; then
 #

diff --git a/tximport_R.R b/tximport_R.R
@@ -2,6 +2,7 @@
 
 library(tximport)
 library(readr)
+library(stringr)
 
 # Rscript tximport_R.R gencode.vM19.metadata.MGI.gz Illumina_PE_SRR.csv output.tsv
 
@@ -10,10 +11,24 @@ args2 = commandArgs(trailingOnly=TRUE)[2]
 args3 = commandArgs(trailingOnly=TRUE)[3]
 
 tx2knownGene <- read_delim(args1, '\t', col_names = c('TXNAME', 'GENEID'))
-exp.table <- read.csv(args2)
+exp.table <- read.csv(args2, row.names=NULL)
 
-files <- paste(c("salmon_output_") , exp.table[,2], c("/quant.sf"), sep='')
-names(files) <- exp.table$name
+files.raw <- exp.table[,2]
+
+# files.raw <- c("SE/test/ttt30.fq.gz", "SE/test/ttt2.fq.gz")
+
+files.raw <- gsub(".gz$", "", files.raw)
+files.raw <- gsub(".fastq$", "", files.raw)
+files.raw <- gsub(".fq$", "", files.raw)
+
+split.vec <- sapply(files.raw, basename)
+# print(paste(c("salmon_output_") , split.vec, c("/quant.sf"), sep=''))
+
+# files <- paste(c("salmon_output_") , exp.table[,2], c("/quant.sf"), sep='')
+files <- paste(c("salmon_output_") , split.vec, c("/quant.sf"), sep='')
+
+
+#print(files)
 
 # txi.salmon <- tximport(files, type = "salmon", tx2gene = tx2knownGene)
 txi.salmon <- tximport(files, type = "salmon", tx2gene = tx2knownGene, countsFromAbundance="scaledTPM")