Last updated: 2018-06-27
workflowr checks: (Click a bullet for more information) ✔ R Markdown file: up-to-date
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
✔ Environment: empty
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
✔ Seed:
set.seed(12345)
The command set.seed(12345)
was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
✔ Repository version: 15c0dd0
wflow_publish
or wflow_git_commit
). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: output/.DS_Store
Untracked files:
Untracked: data/18486.genecov.txt
Untracked: data/YL-SP-18486-T_S9_R1_001-genecov.txt
Untracked: data/bin200.5.T.nuccov.bed
Untracked: data/bin200.Anuccov.bed
Untracked: data/bin200.nuccov.bed
Untracked: data/gene_cov/
Untracked: data/leafcutter/
Untracked: data/nuc6up/
Untracked: data/reads_mapped_three_prime_seq.csv
Untracked: data/ssFC200.cov.bed
Untracked: output/picard/
Untracked: output/plots/
Untracked: output/qual.fig2.pdf
Unstaged changes:
Modified: analysis/dif.iso.usage.leafcutter.Rmd
Modified: analysis/explore.filters.Rmd
Modified: code/Snakefile
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
File | Version | Author | Date | Message |
---|---|---|---|---|
Rmd | 15c0dd0 | Briana Mittleman | 2018-06-27 | try no shift |
html | 8120471 | Briana Mittleman | 2018-06-27 | Build site. |
Rmd | 5e81c8e | Briana Mittleman | 2018-06-27 | add code for macs2 peak calling |
html | ee777df | Briana Mittleman | 2018-06-26 | Build site. |
Rmd | 789d8ef | Briana Mittleman | 2018-06-26 | start test macs analysis. download package |
In this analysis I want to test macs2 as a potential peak caller in the 3’ seq data. This is a widely used peak caller for chip seq data.
I have to create a specific environment to install macs2 because you need to use python 2.7. I call it macs-env. To access this environment I use source activate macs-env.
First, I will merge all of my files into 1 bam file. Using samtools merge. I will do all of my work for this in data/macs2
#!/bin/bash
#SBATCH --job-name=merge
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --partition=broadwl
#SBATCH --mem=20G
#SBATCH --mail-type=END
module load samtools
samtools merge macs2/allBamFiles.bam bam/*.bam
I will create a script in the code directory to call the peaks:
#!/bin/bash
#SBATCH --job-name=macs2nomod
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --output=macs2nomod.out
#SBATCH --error=macs2nomod.err
#SBATCH --partition=broadwl
#SBATCH --mem=8G
#SBATCH --mail-type=END
module load Anaconda3
source activate macs-env
macs2 callpeak -t /project2/gilad/briana/threeprimeseq/data/macs2/allBamFiles.bam -f "BAM" -g 'hs' --outdir /project2/gilad/briana/threeprimeseq/data/macs2/ -n threeprimebatch1 --slocal 1000 --nomodel
This method called 102988 peaks. This is likely more than the true PAS.
Update the -m (MFOLD) term to change the fold enrichment. They must be lower than the upper limmit and higher than the lower limit. The default is 5 50. I will try to make this 20 100.
#!/bin/bash
#SBATCH --job-name=macs2nomod20.100
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --output=macs2nomod.20.100.out
#SBATCH --error=macs2nomod.20.100.err
#SBATCH --partition=broadwl
#SBATCH --mem=8G
#SBATCH --mail-type=END
module load Anaconda3
source activate macs-env
macs2 callpeak -t /project2/gilad/briana/threeprimeseq/data/macs2/allBamFiles.bam -f "BAM" -g 'hs' --outdir /project2/gilad/briana/threeprimeseq/data/macs2/ -n threeprimebatch1.20.100 --slocal 1000 --nomodel -m 20 100
This did not change anything. I am going to try a higher cutoff.
#!/bin/bash
#SBATCH --job-name=macs2nomod40.400
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --output=macs2nomod.40.200.out
#SBATCH --error=macs2nomod.40.200.err
#SBATCH --partition=broadwl
#SBATCH --mem=8G
#SBATCH --mail-type=END
module load Anaconda3
source activate macs-env
macs2 callpeak -t /project2/gilad/briana/threeprimeseq/data/macs2/allBamFiles.bam -f "BAM" -g 'hs' --outdir /project2/gilad/briana/threeprimeseq/data/macs2/ -n threeprimebatch1.40.200 --slocal 1000 --nomodel -m 40 200
Try to not have a shift.
#!/bin/bash
#SBATCH --job-name=macs2noshift
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --output=macs2nomod.noshift.out
#SBATCH --error=macs2nomod.noshift.err
#SBATCH --partition=broadwl
#SBATCH --mem=8G
#SBATCH --mail-type=END
module load Anaconda3
source activate macs-env
macs2 callpeak -t /project2/gilad/briana/threeprimeseq/data/macs2/allBamFiles.bam -f "BAM" -g 'hs' --outdir /project2/gilad/briana/threeprimeseq/data/macs2/ -n threeprimebatch1.noshift --slocal 1000 --nomodel --shift 0 --extsize 200
Does not look like this changed anything.
I also want to run this using seperate files for the total and nuclear fractions.
I wil first merge the total and nuclear bam files seperatly.
#!/bin/bash
#SBATCH --job-name=mergeTN
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --partition=broadwl
#SBATCH --mem=20G
#SBATCH --mail-type=END
module load samtools
samtools merge macs2/TotalBamFiles.bam bam/*T*.bam
samtools merge macs2/NuclearBamFiles.bam bam/*N*.bam
Now I can run the original call peaks on each seperatly with macs2_nomod_TN.sh.
#!/bin/bash
#SBATCH --job-name=macs2nomodTN
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --output=macs2nomodTN.out
#SBATCH --error=macs2nomodTN.err
#SBATCH --partition=broadwl
#SBATCH --mem=8G
#SBATCH --mail-type=END
module load Anaconda3
source activate macs-env
macs2 callpeak -t /project2/gilad/briana/threeprimeseq/data/macs2/TotalBamFiles.bam -f "BAM" -g 'hs' --outdir /project2/gilad/briana/threeprimeseq/data/macs2/ -n Totalthreeprimebatch1 --slocal 1000 --nomodel
macs2 callpeak -t /project2/gilad/briana/threeprimeseq/data/macs2/NuclearBamFiles.bam -f "BAM" -g 'hs' --outdir /project2/gilad/briana/threeprimeseq/data/macs2/ -n Nuclearthreeprimebatch1 --slocal 1000 --nomodel
sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] workflowr_1.0.1 Rcpp_0.12.17 digest_0.6.15
[4] rprojroot_1.3-2 R.methodsS3_1.7.1 backports_1.1.2
[7] git2r_0.21.0 magrittr_1.5 evaluate_0.10.1
[10] stringi_1.2.2 whisker_0.3-2 R.oo_1.22.0
[13] R.utils_2.6.0 rmarkdown_1.8.5 tools_3.4.2
[16] stringr_1.3.1 yaml_2.1.19 compiler_3.4.2
[19] htmltools_0.3.6 knitr_1.18
This reproducible R Markdown analysis was created with workflowr 1.0.1