major change in many scripts

README.md

@@ -19,7 +19,7 @@ RNA-Seq Standard Pipeline is designed to perform standard RNA-Seq analysis for I
 
 4. This pipeline can also perform automatic download from SRA and process the samples. Please see the [template](meta_data/sra_info.txt). Sample names need to be defined in this file for corresponding SRX numberi. The path of [sra_info.txt](meta_data/sra_info.txt) need to be put into [params](params) file. 
 
-5. In [group information file](meta_data/group_info.txt) file, you need to categorize your sample into different groups to perform differential expression analysis.
+5. In [group information file](meta_data/group_info.txt) file, you need to categorize your sample into different groups to perform differential expression analysis. You can have multiple [group information file](meta_data/group_info.txt) with different names.
 
 6. This pipeline use STAR aligner to align the reads. An transcripts annotation GTF format file ("--sjdbGTFfile" option for STAR aligner in the [params](params) file) is required for STAR to extract splice junctions and use them to greatly improve the accuracy of the mapping. And this step also counting number of reads per gene for all the genes in the transcripts annotation GTF file using in the alignment step. The downstream RNA-Seq standard report are all based on this transcript annotation GTF. 
 
@@ -36,9 +36,9 @@ Once everything has been set up, you can run the pipeline in two stages:
 	cut -f1 meta_data/group_info.txt | code/skipn 1 | xargs -n1 -I {} qsub -hard -pe threaded 8 -l tmp_free=190G -l tmp_token=24G -j Y -b Y -cwd -N {} -o {}.out code/By_Sample params {}
 	```
 
-2. After all the samples are processed, you can run the following command to summarize your results:   
+2. After all the samples are processed, you can run the following command to summarize your results (You need to run this command for all your group_info.txt files.):   
 
 	```
-        head -1 meta_data/group_info.txt | cut -f2- | xargs -n1 code/Summarize params pipeline/summarize
+        head -1 meta_data/group_info.txt | cut -f2- | xargs -n1 code/Summarize params pipeline/summarize params
 	```
 3. The summary html report will be generated in the program directory and you need to unzip and open the html file for the report. 

code/By_Sample

@@ -8,12 +8,12 @@ if [ -z $NSLOTS ]
 then
         NSLOTS=6
 fi
-
-SAMPLE=$1
+params=$1
+SAMPLE=$2
 
 FLAG=0
-#if [ -z "$tmp_free" ] || [ -z "$tmp_token" ]
-#then
+if [ ! -z `qstat -j $JOB_ID | grep "tmp_free" | sed 's/ +//g'` ]
+then
 	for CONTENT in `ls -d * | grep -v pipeline`
 	do
 		ln -s $CWD/$CONTENT $TMP/$CONTENT
@@ -30,9 +30,9 @@ FLAG=0
 	done
 	cd $TMP
 	FLAG=1
-#else 
-#	TMP=$CWD
-#fi
+else 
+	TMP=$CWD
+fi
 
 
 
@@ -42,25 +42,24 @@ date
 echo
 
 echo "###### Downloading: $SAMPLE"
-run_fastq-download params pipeline/download meta_data $SAMPLE
-status=$?
+run_fastq-download $params pipeline/download meta_data $SAMPLE
 if [ $FLAG -eq 1 ]; then rsync -rl pipeline/download/ $CWD/pipeline/download/; fi
+status=$?
 if [ $status -ne 0 ]; then exit 1; fi
 
 echo "###### Linking: $SAMPLE"
-run_fastq-link params pipeline/fastq pipeline/download $SAMPLE
+run_fastq-link $params pipeline/fastq pipeline/download $SAMPLE
 status=$?
 if [ $FLAG -eq 1 ]; then rsync -rl pipeline/fastq/ $CWD/pipeline/fastq/; fi
 if [ $status -ne 0 ]; then exit 1; fi
-
 echo "###### Alignning: $SAMPLE"
-run_alignment params pipeline/alignment pipeline/fastq $SAMPLE
+run_alignment $params pipeline/alignment pipeline/fastq $SAMPLE
 status=$?
 if [ $FLAG -eq 1 ]; then rsync -rl pipeline/alignment/ $CWD/pipeline/alignment/; fi
 if [ $status -ne 0 ]; then exit 1; fi
 
 echo "###### Assemblying: $SAMPLE"
-run_assembly params pipeline/cufflinks pipeline/alignment $SAMPLE
+run_assembly $params pipeline/cufflinks pipeline/alignment $SAMPLE
 status=$?
 if [ $FLAG -eq 1 ]; then rsync -rl pipeline/cufflinks/ $CWD/pipeline/cufflinks/; fi
 if [ $status -ne 0 ]; then exit 1; fi

code/DESeq2.r

@@ -6,7 +6,7 @@ Sys.time()
 args <- commandArgs(trailingOnly = TRUE)
 if (length(args)!=5) { write(usage,stderr()); quit(save='no'); }
 
-my.packages=c("DESeq2","GenomicFeatures","dplyr","BiocParallel")
+my.packages=c("DESeq2","GenomicFeatures","dplyr","BiocParallel","tools")
 for (p in my.packages){
     suppressPackageStartupMessages(library(p,character.only=TRUE,verbose=FALSE))
 }
@@ -20,6 +20,7 @@ id2name=as.matrix(read.table(args[2],header=F,row.names = 1))
 id2name=as.matrix(id2name[order(rownames(id2name)),])
 group=read.table(args[3],header=T,fill=T)
 group_by=args[4]
+outdir=file_path_sans_ext(basename(args[3]))
 grp=factor(group[,group_by])
 countsData <-read.table(args[5],header=T,row.names=1)
 # raw count must be sorted by rowname before run DESeq2 because when matching with GTF file, it doesn't sort automaticlly and will cause mistake
@@ -37,43 +38,43 @@ dds <- DESeqDataSetFromMatrix(countsData, colData=data.frame(sampleData), design
 # add GTF file into DESeq2 Matrix
 rowRanges(dds) <- ebg
 dds <- DESeq(dds)
-
 # Regularized-logarithm transformation
 rld <- rlog(dds, blind=FALSE)
 rownames(rld)=cbind(id2name[,1][match(rownames(rld),rownames(id2name))],rld)[,1]
 
+dir.create(file.path(outdir,group_by),recursive=T,showWarnings=F)
+
 # extract norm count
 norm_count=counts(dds, normalized=T)
 # add gene name
 norm_count=data.frame(cbind(Gene=id2name[,1][match(rownames(norm_count),rownames(id2name))],norm_count))
 # convert to Numeric for the count field, for ReportingTools
-if(!file.exists("norm_count.txt")){write.table(norm_count,file="norm_count.txt",quote=F,sep="\t",row.names=F,col.names=T)}
+if(!file.exists(file.path(outdir,"norm_count.txt"))){write.table(norm_count,file=file.path(outdir,"norm_count.txt"),quote=F,sep="\t",row.names=F,col.names=T)}
 rownames(norm_count)=norm_count$Gene
 norm_count=norm_count[,-1]
 
+
 # extract fpkm 
 fpkm_table=fpkm(dds,robust=T)
 # add gene name
 fpkm_table=data.frame(cbind(Gene=id2name[,1][match(rownames(fpkm_table),rownames(id2name))],fpkm_table))
 # convert to Numeric for the fpkm field, for ReportingTools
-if(!file.exists("fpkm_table.txt")){write.table(fpkm_table,file="fpkm_table.txt",quote=F,sep="\t",row.names=F,col.names=T)}
+if(!file.exists(file.path(outdir,"fpkm_table.txt"))){write.table(fpkm_table,file=file.path(outdir,"fpkm_table.txt"),quote=F,sep="\t",row.names=F,col.names=T)}
 rownames(fpkm_table)=fpkm_table$Gene
 fpkm_table=fpkm_table[,-1]
 
 comparisons=cbind("grp",t(combn(levels(colData(dds)$grp),2)))
 
 ToCharacter=function(x){x}
 
-dir.create(group_by)
-
 for (i in 1:dim(comparisons)[1]){
   res=results(dds,contrast=comparisons[i,])
   id2name=as.matrix(id2name[order(rownames(id2name)),])
   res$symbol=apply(id2name,1,ToCharacter)
   mytable <- res %>% data.frame() %>% select(symbol,log2FoldChange,padj,pvalue) %>% mutate(log2FoldChange=round(log2FoldChange,digits=4),padj=round(padj,digits=4),pvalue=round(pvalue,digits=4))
-  colnames(mytable)=c("_ID_",paste(comparisons[i,2],"_vs_",comparisons[i,3],"_FC",sep=""),paste(comparisons[i,2],"_vs_",comparisons[i,3],"_FDR",sep=""),paste(comparisons[i,2],"_vs_",comparisons[i,3],"_PVAL",sep=""))
-  write.table(mytable,file=paste(group_by,"/",comparisons[i,2],"_vs_",comparisons[i,3],"_DE.txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
+  colnames(mytable)=c("_ID_",paste(comparisons[i,2],"_vs_",comparisons[i,3],"_LFC",sep=""),paste(comparisons[i,2],"_vs_",comparisons[i,3],"_FDR",sep=""),paste(comparisons[i,2],"_vs_",comparisons[i,3],"_PVAL",sep=""))
+  write.table(mytable,file=paste(outdir,"/",group_by,"/",comparisons[i,2],"_vs_",comparisons[i,3],"_DE.txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
 }
 
-save.image(paste(group_by,"/",group_by,"_Rscript.RData",sep=""))
+save.image(paste0(outdir,"/",group_by,"/",group_by,"_Rscript.RData"))
 Sys.time()

code/DESeq_post_v3.r

@@ -23,9 +23,11 @@ toNumeric<-function(myMatrix){
 fpkm_table=toNumeric(fpkm_table)
 norm_count=toNumeric(norm_count)
 
-outdir=basename(dirname(dataset))
-dir.create(outdir)
-htmlRep <- HTMLReport(shortName = outdir, title="RNA-Seq Experiment Standard Report", reportDirectory = outdir)
+longdir=dirname(dataset)
+
+if(grep("/",longdir)){outdir=basename(longdir)}else{outdir=longdir}
+dir.create(longdir,showWarnings=F)
+htmlRep <- HTMLReport(shortName = outdir, title="RNA-Seq Experiment Standard Report", reportDirectory = longdir)
 
 
 # signature as goldern standard to be tested
@@ -53,7 +55,7 @@ align_reads <- function(x){
     scale_fill_manual(values=c("limegreen","forestgreen","orchid","indianred1")) +
     ggtitle("b. Alignment Category Reads") +
     xlab("") + ylab("Number of Reads in Millions") +
-    theme(axis.title = element_text(size=8),axis.text.y = element_blank(),axis.text.x = element_text(size=5),plot.title = element_text(size=8),legend.title=element_blank(),legend.text=element_text(size=8),axis.ticks.y=element_blank(), axis.line.y=element_blank()) +
+    theme(axis.title = element_text(size=8),axis.text.y = element_blank(),plot.title = element_text(size=8),legend.title=element_blank(),legend.text=element_text(size=8),axis.ticks.y=element_blank(), axis.line.y=element_blank()) +
     coord_flip()
     return(myplot)
 }
@@ -68,7 +70,7 @@ align_PCT <- function(x){
     ggtitle("a. Alignment Category Percentage") +
     scale_fill_manual(values=c("limegreen","forestgreen","orchid","indianred1")) +
     xlab("") + ylab("Percentage of Reads") +
-    theme(axis.title = element_text(size=8),axis.text.y = element_text(size=2.8,angle=0),axis.text.x = element_text(size=5),plot.title = element_text(size=8)) +
+    theme(axis.title = element_text(size=8),axis.text.y = element_text(angle=0),plot.title = element_text(size=8)) +
     guides(fill=FALSE) +
     coord_flip()
     return(myplot)
@@ -92,10 +94,10 @@ splDistPlot <- function(mymatrix){
   rownames(sampleDistMatrix) <- colnames(mymatrix)
   colnames(sampleDistMatrix) <- colnames(mymatrix)
   colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
-  png(filename=paste(outdir,"/figures",outdir,"/Sample_Distance.png",sep=""))
+  png(filename=paste(longdir,"/figures",outdir,"/Sample_Distance.png",sep=""))
   pheatmap(sampleDistMatrix, clustering_distance_rows=sampleDists, clustering_distance_cols=sampleDists, col=colors, main="Distance Between Samples")
   dev.off()
-  pdf(file=paste(outdir,"/figures",outdir,"/Sample_Distance.pdf",sep=""),onefile=F)
+  pdf(file=paste(longdir,"/figures",outdir,"/Sample_Distance.pdf",sep=""),onefile=F)
   pheatmap(sampleDistMatrix, clustering_distance_rows=sampleDists, clustering_distance_cols=sampleDists, col=colors, main="Distance Between Samples")
   dev.off()
 }
@@ -146,27 +148,27 @@ geneBarplot <- function(mytitle,data){
 
 geneBoxplot <- function(mytitle,data){
   if(dim(data)[1]!=0){
-    myplot <- ggplot(mydata, aes(x=grp, y=count)) + 
-    geom_boxplot(aes(fill=grp)) + 
-    geom_jitter(size=0.1,width=0.1) +
-    facet_grid(gene~.) +
-    theme(strip.text.y = element_text(size=8, angle=0), legend.title=element_blank(), axis.text.y = element_blank() , axis.ticks.y = element_blank()) +
-    labs(title = paste("FPKM for ", mytitle," Signature Genes",sep="")) +
-    guides(color=FALSE) + 
-    xlab("Gene") + ylab("FPKM") +
-    coord_flip()
+    myplot <- ggplot(data, aes(x=grp, y=count)) +
+      geom_boxplot(aes(fill=grp,color=grp)) +
+      facet_grid(gene~.) +
+      theme(strip.text.y = element_text(size=8, angle=0), legend.title=element_blank(), axis.text.y = element_blank() , axis.ticks.y = element_blank(), axis.line.y = element_blank()) +
+      labs(title = paste("FPKM for ", mytitle," Signature Genes",sep="")) +
+      guides(color=FALSE) +
+      xlab("Gene") + ylab("FPKM") +
+      coord_flip()+ 
+      stat_summary(geom = "crossbar", fatten=0, color="white", fun.data = function(x){ return(c(y=median(x), ymin=median(x), ymax=median(x))) })
     return(myplot)
   }else{return(NULL)}
 }
 
-highexp <- mydata %>% group_by(gene) %>% filter(max(count)>100)
+highexp <- mydata %>% group_by(gene) %>% filter(max(count)>100) %>% ungroup
 mytitle="High Exp"
-myplot <- geneBarplot(mytitle,highexp)
+myplot <- geneBoxplot(mytitle,highexp)
 if(!is.null(myplot)){publish(myplot, htmlRep,name=paste("FPKM for ",mytitle," Signature Genes",sep=""))}
 
-lowexp <- mydata %>% group_by(gene) %>% filter(max(count)<100)
+lowexp <- mydata %>% group_by(gene) %>% filter(max(count)<100) %>% ungroup
 mytitle="Low Exp"
-myplot <- geneBarplot(mytitle,lowexp)
+myplot <- geneBoxplot(mytitle,lowexp)
 if(!is.null(myplot)){publish(myplot, htmlRep,name=paste("FPKM for ",mytitle," Signature Genes",sep=""))}
 
 publish("FPKM Table:", htmlRep)
@@ -192,10 +194,10 @@ publish(hwrite("Download", link = paste("../raw_count.txt",sep="")), htmlRep)
 #plot(res$log2FoldChange,-log(res$padj,10),main="Volcano Plot of FDR-adjusted P-values",pch=20,cex=0.4,xlab="log2(FC of mean of normalized counts)",ylab="-log10(adj P-value)")
 
 # Dispersion Plot. First, gene-wise MLEs are obtained using only the respective gene’s data (black dots). Then, a curve (red) is fit to the MLEs to capture the overall trend of dispersion-mean dependence. This fit is used as a prior mean for a second estimation round, which results in the final MAP estimates of dispersion (arrow heads). This can be understood as a shrinkage (along the blue arrows) of the noisy gene-wise estimates toward the consensus represented by the red line. The black points circled in blue are detected as dispersion outliers and not shrunk toward the prior (shrinkage would follow the dotted line). 
-png(filename=paste(outdir,"/figures",outdir,"/Dispersion_Plot.png",sep=""))
+png(filename=paste(longdir,"/figures",outdir,"/Dispersion_Plot.png",sep=""))
 plotDispEsts(dds, main="Dispersion Plot")
 dev.off()
-pdf(file=paste(outdir,"/figures",outdir,"/Dispersion_Plot.pdf",sep=""))
+pdf(file=paste(longdir,"/figures",outdir,"/Dispersion_Plot.pdf",sep=""))
 plotDispEsts(dds, main="Dispersion Plot")
 dev.off()
 himg <- hwriteImage(paste("figures",outdir,"/Dispersion_Plot.png",sep=""),link=paste("figures",outdir,"/Dispersion_Plot.pdf",sep=""))
@@ -206,7 +208,7 @@ publish(paste("Analysis for Different Comparison Groups ",sep=""), htmlRep)
 for (i in 1:dim(comparisons)[1]){
   preprocessing="rld"
   print(comparisons[i,])
-  myreport=summarize_comparison(outdir,signature,target,"DE",comparisons[i,],preprocessing,fpkm_table)
+  myreport=summarize_comparison(longdir,signature,target,"DE",comparisons[i,],preprocessing,fpkm_table)
   publish(hwrite(myreport, link = paste(myreport,"/","report.html",sep="")), htmlRep)
 }
 

code/Summarize

@@ -14,16 +14,25 @@ elif [[ $NSLOTS == "" ]]; then NSLOTS=6; fi
 
 PARAMS=$1
 OUT_DIR=$2
-GROUP=$3
+group_info=$3
+GROUP=$4
+group_sheet=`basename $group_info | sed 's/\.txt$//g'`
 
 params=`cat $PARAMS | grep "^DESeq2" | cut -f2`
 
 cd $OUT_DIR
 
+mkdir -p $group_sheet
 
-align_summary alignment/*/STAR.Log.final.out > align_summary.txt
+if [ ! -f $group_sheet/align_summary.txt ]
+then
+align_summary `cat $group_info | cut -f1 | skipn 1 | sed 's/^/alignment\//g' | sed 's/$/\/STAR.Log.final.out/g' | tr '\n' ' '` > $group_sheet/align_summary.txt
+fi
 
-summarize_raw_count alignment/*/STAR.count.txt > raw_count.txt
+if [ ! -f $group_sheet/raw_count.txt ]
+then
+summarize_raw_count `cat $group_info | cut -f1 | skipn 1 | sed 's/^/alignment\//g' | sed 's/$/\/STAR.count.txt/g' | tr '\n' ' '` > $group_sheet/raw_count.txt
+fi
 
 #qsub -b Y -cwd -j Y -pe threaded 8 /local/apps/R/3.2.0/bin/Rscript ~/utilities/yixiao_scripts/count2fpkm.r ../referenceFiles/Homo_sapiens.GRCh37.\?\?.gtf ../referenceFiles/Homo_sapiens.GRCh37.\?\?.id2name.txt ../alignment/raw_count.txt
 
@@ -35,12 +44,12 @@ then
 	gtf_id2name ${params}
 fi 
 
-DESeq2.r $params $ID2NAME meta_data/group_info.txt $GROUP raw_count.txt
+DESeq2.r $params $ID2NAME $group_info $GROUP $group_sheet/raw_count.txt
 
-DESeq_post_v3.r $GROUP/${GROUP}_Rscript.RData align_summary.txt meta_data/signature.txt meta_data/target.txt
+DESeq_post_v3.r $group_sheet/$GROUP/${GROUP}_Rscript.RData $group_sheet/align_summary.txt meta_data/signature.txt meta_data/target.txt
 
-rm -rf align_summary.txt Rplots.pdf
+rm -rf Rplots.pdf
 
-zip -r ${GROUP}_Report.zip $GROUP fpkm_table.txt norm_count.txt raw_count.txt
+#zip -r ${group_sheet}_Report.zip $group_sheet
 
-ln -s $OUT_DIR/${GROUP}_Report.zip ../../${GROUP}_Report.zip
+#ln -s $OUT_DIR/${group_sheet}_Report.zip ../../${group_sheet_Report}.zip

code/install_R_packages.r

@@ -1,6 +1,6 @@
 #!/usr/bin/env Rscript
 
-list.of.packages=c("DESeq2","GenomicFeatures","dplyr","BiocParallel","pheatmap","RColorBrewer","ggplot2","ReportingTools","hwriter","reshape2","preprocessCore","cowplot","diagram")
+list.of.packages=c("tools","DESeq2","GenomicFeatures","dplyr","BiocParallel","pheatmap","RColorBrewer","ggplot2","ReportingTools","hwriter","reshape2","preprocessCore","cowplot","diagram")
 missing.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
 for (p in missing.packages){
  if(!file.exists(Sys.getenv("R_LIBS_USER"))){

code/run_assembly

@@ -64,7 +64,11 @@ cufflinks -p $NSLOTS --library-type=fr-${STRAND} $params -o $BASE_DIR $INP_DIR/$
 "
 #printf "\n\n $CMD \n\n"
 printf "\n\n $CMD \n\n" >> $BASE_DIR/run.bash
-eval $CMD
+
+
+source code/custom_bash_profile
+cufflinks -p $NSLOTS --library-type=fr-${STRAND} $params -o $BASE_DIR $INP_DIR/$SAMPLE/accepted_hits.bam
+
 
 fi
 

code/run_fastq-link

@@ -22,7 +22,7 @@ date
 #########################
 # input
 PARAMS=$1
-OUT_DIR=$2
+OUT_DIR=`echo $2 | sed 's/\/\//\//g'`
 INP_DIR=$3
 SAMPLE=$4
 
@@ -42,19 +42,20 @@ then
 for FILE in $params
 do
 	mylink=`head -1 $FILE`
+echo $mylink
 	if [ -d $mylink ]
 	then
 		dirname=`basename $FILE | sed 's/.txt//g'`
 		ln -s $mylink $OUT_DIR/$dirname
 		for record in `cat $FILE | grep -v "\/"`
 		do
 			sample_name=`echo $record | cut -f1`		
-			if [ ! -d $OUT_DIR/$sample_name ]
-			then
-				mkdir -p $OUT_DIR/$sample_name
-			fi
 			if [ "$sample_name" == "$SAMPLE" ] || [ -z "$SAMPLE" ]
 			then
+				if [ ! -d $OUT_DIR/$sample_name ]
+				then
+					mkdir -p $OUT_DIR/$sample_name
+				fi
 				identifier=`echo $record | cut -f2`
 				printf "\n\n***** Linking $sample_name To $mylink/${identifier}*...\n"
 				for myfastq in `ls $OUT_DIR/$dirname/* | grep "$identifier"`

code/summarize_raw_count

@@ -10,4 +10,4 @@ do
 	INDEX=${INDEX}","${i}
 done
 
-paste $@ | cut -f${INDEX} > raw_count.txt
+paste $@ | cut -f${INDEX} 

params

@@ -1,9 +1,9 @@
 #STAR-Genome	--genomeFastaFiles referenceFiles/genome.fa 
 fastq-download	meta_data
-fastq-link	meta_data/20160224.txt #meta_data/sra_info.txt
+fastq-link	meta_data/20160425.txt #meta_data/sra_info.txt
 #fastqc	meta_data/sra_info.txt
 STAR	--readFilesCommand zcat --outReadsUnmapped Fastx --outSAMtype BAM SortedByCoordinate --outFilterMultimapNmax 1 --genomeDir referenceFiles/STAR_Reference --sjdbGTFfile referenceFiles/gencode.v19.annotation.gtf
 MarkDuplicates	ASSUME_SORTED=false CREATE_INDEX=true REMOVE_DUPLICATES=true
 rnaseq_tracks	referenceFiles/chromInfo.txt
-#cufflinks	-q -M referenceFiles/maskPlussnRNA.gtf -g referenceFiles/gencode.v19.annotation.gtf
+cufflinks	-q -M referenceFiles/maskPlussnRNA.gtf -g referenceFiles/gencode.v19.annotation.gtf
 DESeq2	referenceFiles/gencode.v19.annotation.gtf