Dear Bo,

Sounds like you have not clustered reads into peaks yet. After you do, you should not have too many distinct potential polyA sites. FYI, we got 205470 peaks for one time point and 75511 for another time point. Sarah Sheppard, a Ph.D student who recently graduated, wrote some perl script to do such clustering. I cced her in this email in case you are interested in her scripts.

Also, we divide all potential peaks into batches, each with about 10,000 peaks and run them in the school cluster hpcc  simultaneously.

Here is how to run your R script using qsub in the cluster assuming your R script is called predPas.R .

qsub –cwd predPas.sh

 Here is the content in the predPas.sh

#!/bin/bash
#$ -V
#$ -l mem_free=26G
/share/bin/R-3.0.1/lib64/R/bin/R CMD BATCH predPas.R

The R session output will be in predPas.Rout. For missing libraries, you could contact the cluster management group to install the proper R packages.

Here is the user manual for cleanUpdTSeq if you have not read it http://www.bioconductor.org/packages/2.13/bioc/vignettes/cleanUpdTSeq/inst/doc/cleanUpdTSeq.pdf. We probably need to add a section on how to run it in a cluster environment and add the clustering script.

Thanks for your feedback!

Best regards,

Julie



On 9/12/13 10:04 AM, "Han, Bo" <Bo.Han@umassmed.edu> wrote:

Hello, Dear Prof. Zhu,

This is a Graduate Student in the Zamore Lab.

Congratulations on your recent paper on bioinformatics and thanks to that, our analysis on PAS could be easier.

Xin, my collaborator has met with you last week and got the document of cleanUpdTSeq. I was following that document, but met a problem.

Seems that the program requires lots of memory and I couldn't finish the running even on our best node in school. I have merged the reads that mapped to the same position, ,and even separated the reads to different chromosomes (each 5—10M). But still took too much memory and cannot finish running.

I would really appreciate it if you can give me some suggestions. The data (with individual chromosomes) are located at /home/hanb/scratch/prepachytene/PAS/mouse.testis.wild.type.adult.PASseq/chromosomes.

Here is the code that I ran the R:

#source("http://bioconductor.org/biocLite.R")
#biocLite("cleanUpdTSeq")
library(cleanUpdTSeq)
library(BSgenome.Mmusculus.UCSC.mm9)

args = commandArgs (T)
inputBedFile = args[1]

testSet = read.table(inputBedFile, sep = "\t", header = F)

peaks = BED2GRangesSeq (testSet, withSeq=F)

testSet.NaiveBayes = buildFeatureVector (peaks,BSgenomeName = Mmusculus, upstream = 40,
downstream = 30, wordSize = 6, alphabet=c("ACGT"),
sampleType = "unknown",replaceNAdistance = 30,
method = "NaiveBayes", ZeroBasedIndex = 1, fetchSeq = TRUE)

data(data.NaiveBayes)

predictTestSet(data.NaiveBayes$Negative, data.NaiveBayes$Positive, testSet.NaiveBayes,
outputFile = paste (inputBedFile, ".tsv", sep=""), assignmentCutoff = 0.5)

And how I got the bed from fastq and ran the Rscript (last two lines)

#! /bin/bash -x

################
# Major Config #
################
export PPP_PATH=/home/hanb/nearline/ppp
export PATH=$PPP_PATH:$PATH

# FQ input
FQ=$1
# CPU
CPU=8
# step
STEP=1
# genome
Genome=mm9
# bowtie2 index
BOWTIE2_INDEX=/home/hanb/nearline/Mus_musculus/UCSC/mm9/Sequence/WholeGenomeFasta/genome
# chrom
Chrome=/home/hanb/nearline/small_RNA_Pipeline/common_files/mm9/ChromInfo.txt
# GTF file with gene annotation
GTF_FILE=/home/hanb/scratch/mm9.newest.07dpp.gtf

# check read length
READ_LEN=`head -2 $FQ | awk '{getline; printf "%d", length($1)}'`
# STAR sjdbOverhang option
STARsjdbOverhang=$((READ_LEN-1))
# STAR genome index directory
STAR_Genome_INDEX_DIR=/home/hanb/nearline/PE_Pipeline/common_files/mm9/STARgenomeIndex
STAR_rRNA_INDEX_DIR=/home/hanb/nearline/PE_Pipeline/common_files/mm9/rRNAIndex
# STAR genome index
STAR_Genome_INDEX=$STAR_Genome_INDEX_DIR/$STARsjdbOverhang
STAR_rRNA_INDEX=$STAR_rRNA_INDEX_DIR/$STARsjdbOverhang
# MM
MM=$(((READ_LEN+1)/25))

# prefix
PREFIX=${FQ%.fq}
# trimmed 3' adaptor and filter quality
ADAPTOR_SEQ="AATGGAATTC" # allow two As in front of the 3' adaptor
# ADAPTOR_SEQ="AAAAAAAAAAAAAAAAAAAATGGAATTCTCGGGTGCCAAGGC"
# reads are subjected to quality trimming first until reaching phred >= 3
# then find and trim adaptor
# reads without adaptor found are still reported but later filtered
# reads with N is dumped
# reads shorter than 25nt after quality trimming and adaptor clipping are dumped as well
[ ! -f .status.${STEP}.flexbar ] && \
flexbar \
-r $FQ \
-t ${PREFIX}.trimmed \
-f fastq-i1.8 \
-n $CPU \
-as $ADAPTOR_SEQ \
-at 1.0 \
-ao 10 \
--min-readlength 25 \
--max-uncalled 0 \
-q 3 && \
awk '{NAME=$1;getline;SEQ=$1;getline;getline;QUAL=$1; ++A; if (substr (SEQ,length($0)-5)=="AAAAAA") {++B; printf "%s\n%s\n+\n%s\n", NAME, SEQ, QUAL}}END{print A"\t"B >> "filter_6A.log"}' ${PREFIX}.trimmed.fastq > ${PREFIX}.trimmed.fq && \
touch .status.${STEP}.flexbar
STEP=$((STEP+1))

# map to genome using STAR with stringent condition
# ideally, polyA by internal priming shoud be able to be mapped
INPUT=${PREFIX}.trimmed.fq
[ ! -f .status.${STEP}.genome_mapping_with_polyA ] && \
STAR \
--runMode alignReads \
--genomeDir $STAR_Genome_INDEX \
--readFilesIn ${INPUT} \
--runThreadN $CPU \
--outFilterScoreMin 0 \
--outFilterScoreMinOverLread 0.90 \
--outFilterMatchNmin 0 \
--outFilterMatchNminOverLread 0.90 \
--outFilterMultimapScoreRange 1 \
--outFilterMultimapNmax -1 \
--outFilterMismatchNmax 10 \
--outFilterMismatchNoverLmax 0.10 \
--alignIntronMax 0 \
--alignIntronMin 21 \
--alignSJDBoverhangMin 1 \
--outFilterIntronMotifs RemoveNoncanonicalUnannotated \
--outSAMstrandField intronMotif \
--outSAMattributes All \
--genomeLoad NoSharedMemory \
--outFileNamePrefix ${PREFIX}.x_rRNA.${Genome}.+polyA.unmapped.softSlipping. \
--outReadsUnmapped Fastx \
--outSJfilterReads Unique \
--seedSearchStartLmax 20 \
--seedSearchStartLmaxOverLread 1.0 \
--chimSegmentMin 0 && \
touch .status.${STEP}.genome_mapping_with_polyA
STEP=$((STEP+1))

# map unmapped without polyA
## run trimpoly
## convert FQ to FA for trimpoly

[ ! -f .status.${STEP}.trimpoly ] && \
awk '{name=substr($1,2);print ">"name; getline;print $1;getline;getline}' ${PREFIX}.trimmed.fq > ${PREFIX}.trimmed.fa && \
trimpoly -e 3 -s 5 < ${PREFIX}.trimmed.fa > ${PREFIX}.trimmed.trimpoly.log &&
trimpoly2 ${PREFIX}.trimmed.fq ${PREFIX}.trimmed.trimpoly.log > ${PREFIX}.trimmed.-polyA.fq && \
touch .status.${STEP}.trimpoly
STEP=$((STEP+1))

# map to genome again
INPUT=${PREFIX}.trimmed.-polyA.fq
[ ! -f .status.${STEP}.genome_mapping_without_polyA ] && \
STAR \
--runMode alignReads \
--genomeDir $STAR_Genome_INDEX \
--readFilesIn ${INPUT} \
--runThreadN $CPU \
--outFilterScoreMin 0 \
--outFilterScoreMinOverLread 0.72 \
--outFilterMatchNmin 0 \
--outFilterMatchNminOverLread 0.72 \
--outFilterMultimapScoreRange 1 \
--outFilterMultimapNmax -1 \
--outFilterMismatchNmax 15 \
--outFilterMismatchNoverLmax 0.15 \
--alignIntronMax 0 \
--alignIntronMin 21 \
--alignSJDBoverhangMin 1 \
--outFilterIntronMotifs RemoveNoncanonicalUnannotated \
--outSAMstrandField intronMotif \
--outSAMattributes All \
--genomeLoad NoSharedMemory \
--outFileNamePrefix ${PREFIX}.x_rRNA.${Genome}.-polyA. \
--outReadsUnmapped Fastx \
--outSJfilterReads Unique \
--seedSearchStartLmax 20 \
--seedSearchStartLmaxOverLread 1.0 \
--chimSegmentMin 0 && \
touch .status.${STEP}.genome_mapping_without_polyA
STEP=$((STEP+1))

# getting statistics
InputReads=`grep 'Number of input reads' ${PREFIX}.x_rRNA.${Genome}.-polyA.Log.final.out | awk '{print $NF}'`
UniquReads=`grep 'Uniquely mapped reads number' ${PREFIX}.x_rRNA.${Genome}.-polyA.Log.final.out | awk '{print $NF}'`
MultiReads=`grep 'Number of reads mapped to multiple loci' ${PREFIX}.x_rRNA.${Genome}.-polyA.Log.final.out | awk '{print $NF}'`
AllMapReads=$((UniquReads+MultiReads))
UnMapReads=$((InputReads-UniquReads-MultiReads))
# normalization factor
NormScale=`echo $UniquReads | awk '{printf "%f",1000000.0/$1}'`

# convert format
[ ! -f .status.${STEP}.data_convert1 ] && \
samtools view -bS ${PREFIX}.x_rRNA.${Genome}.-polyA.Aligned.out.sam > ${PREFIX}.x_rRNA.${Genome}.-polyA.Aligned.out.bam && \
samtools sort -@ $CPU ${PREFIX}.x_rRNA.${Genome}.-polyA.Aligned.out.bam ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted && \
rm -rf ${PREFIX}.x_rRNA.${Genome}.-polyA.Aligned.out.sam ${PREFIX}.x_rRNA.${Genome}.-polyA.Aligned.out.bam && \
touch .status.${STEP}.data_convert1
STEP=$((STEP+1))

# making bigWig
[ ! -f .status.${STEP}.bigWig ] && \
samtools view -h ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.bam | \
filterSoftClippedReads -3 | \
samtools view -bS - | \
bedtools bamtobed -i - | \
awk 'BEGIN{OFS="\t";FS="\t"}{if ($6=="+") {$2=$3-1} else {$3=$2+1} print $0}' | \
uniq \
> ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.bed && \
awk 'BEGIN{OFS="\t"}{if ($5==255) {print $0 > "/dev/stdout"} else {print $0 > "/dev/stderr"}}' ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.bed \
1> ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.bed \
2> ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.multip.bed && \
paraFile="${RANDOM}${RANDOM}.para" && \
echo "bedtools genomecov -scale $NormScale -bg -strand + -i ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.bed -g $Chrome > ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.plus.bedgraph && bedGraphToBigWig ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.plus.bedgraph $Chrome ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.plus.BW"  > $paraFile && \
echo "bedtools genomecov -scale $NormScale -bg -strand - -i ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.bed -g $Chrome | awk 'BEGIN{OFS=\"\t\"}{\$4 = -\$4; print \$0}' > ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.minus.bedgraph && bedGraphToBigWig ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.minus.bedgraph $Chrome ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.minus.BW" >> $paraFile && \
ParaFly -CPU $CPU -c $paraFile -failed_cmds ${STEP}.failed_commands && \
echo "track name=U.5p.${PREFIX}(+) description=${PREFIX}(+) maxHeightPixels=25 alwaysZero=on autoScale=on  yLineMark=0 yLineOnOff=on type=bigWig color=$((RANDOM%255)),$((RANDOM%255)),$((RANDOM%255)) visibility=full bigDataUrl=http://zlab.umassmed.edu/~hanb/Bo/${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.plus.BW"  >  ${PREFIX}.x_rRNA.${Genome}.-polyA.5p.unique.track && \
echo "track name=U.5p.${PREFIX}(-) description=${PREFIX}(-) maxHeightPixels=25 alwaysZero=on autoScale=on  yLineMark=0 yLineOnOff=on type=bigWig color=$((RANDOM%255)),$((RANDOM%255)),$((RANDOM%255)) visibility=full bigDataUrl=http://zlab.umassmed.edu/~hanb/Bo/${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.minus.BW" >> ${PREFIX}.x_rRNA.${Genome}.-polyA.5p.unique.track && \
touch .status.${STEP}.bigWig
STEP=$((STEP+1))


[ ! -f .status.${STEP}.doCleanUpdTSeq ] && \
awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$2,$3,1,1,$6}' ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.unique.bed | uniq | awk 'BEGIN{FS="\t";OFS="\t"}{print $1,$2,$3,NR,1,$6}' > ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.bed2 && \
Rscript $PPP_PATH/cleanUpdTSeq.R ${PREFIX}.x_rRNA.${Genome}.-polyA.sorted.noS.5p.bed2 && \
touch .status.${STEP}.doCleanUpdTSeq










On Sep 4, 2013, at 7:28 PM, "Li, Xin" <Xin.Li@umassmed.edu> wrote:

<cleanUpdTSeq.pdf>



	[[alternative HTML version deleted]]

