Last updated: 2019-02-15
Checks: 6 0
Knit directory: threeprimeseq/analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.2.0). The Report tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
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.
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.
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.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use 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: data/.DS_Store
Ignored: data/perm_QTL_trans_noMP_5percov/
Ignored: output/.DS_Store
Untracked files:
Untracked: KalistoAbundance18486.txt
Untracked: analysis/4suDataIGV.Rmd
Untracked: analysis/DirectionapaQTL.Rmd
Untracked: analysis/EvaleQTLs.Rmd
Untracked: analysis/YL_QTL_test.Rmd
Untracked: analysis/ncbiRefSeq_sm.sort.mRNA.bed
Untracked: analysis/snake.config.notes.Rmd
Untracked: analysis/verifyBAM.Rmd
Untracked: analysis/verifybam_dubs.Rmd
Untracked: code/PeaksToCoverPerReads.py
Untracked: code/strober_pc_pve_heatmap_func.R
Untracked: data/18486.genecov.txt
Untracked: data/APApeaksYL.total.inbrain.bed
Untracked: data/ApaQTLs/
Untracked: data/ChromHmmOverlap/
Untracked: data/DistTXN2Peak_genelocAnno/
Untracked: data/GM12878.chromHMM.bed
Untracked: data/GM12878.chromHMM.txt
Untracked: data/LianoglouLCL/
Untracked: data/LocusZoom/
Untracked: data/NuclearApaQTLs.txt
Untracked: data/PeakCounts/
Untracked: data/PeakCounts_noMP_5perc/
Untracked: data/PeakCounts_noMP_genelocanno/
Untracked: data/PeakUsage/
Untracked: data/PeakUsage_noMP/
Untracked: data/PeakUsage_noMP_GeneLocAnno/
Untracked: data/PeaksUsed/
Untracked: data/PeaksUsed_noMP_5percCov/
Untracked: data/RNAkalisto/
Untracked: data/RefSeq_annotations/
Untracked: data/TotalApaQTLs.txt
Untracked: data/Totalpeaks_filtered_clean.bed
Untracked: data/UnderstandPeaksQC/
Untracked: data/WASP_STAT/
Untracked: data/YL-SP-18486-T-combined-genecov.txt
Untracked: data/YL-SP-18486-T_S9_R1_001-genecov.txt
Untracked: data/YL_QTL_test/
Untracked: data/apaExamp/
Untracked: data/apaQTL_examp_noMP/
Untracked: data/bedgraph_peaks/
Untracked: data/bin200.5.T.nuccov.bed
Untracked: data/bin200.Anuccov.bed
Untracked: data/bin200.nuccov.bed
Untracked: data/clean_peaks/
Untracked: data/comb_map_stats.csv
Untracked: data/comb_map_stats.xlsx
Untracked: data/comb_map_stats_39ind.csv
Untracked: data/combined_reads_mapped_three_prime_seq.csv
Untracked: data/diff_iso_GeneLocAnno/
Untracked: data/diff_iso_proc/
Untracked: data/diff_iso_trans/
Untracked: data/ensemble_to_genename.txt
Untracked: data/example_gene_peakQuant/
Untracked: data/explainProtVar/
Untracked: data/filtPeakOppstrand_cov_noMP_GeneLocAnno_5perc/
Untracked: data/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.bed
Untracked: data/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.noties.bed
Untracked: data/first50lines_closest.txt
Untracked: data/gencov.test.csv
Untracked: data/gencov.test.txt
Untracked: data/gencov_zero.test.csv
Untracked: data/gencov_zero.test.txt
Untracked: data/gene_cov/
Untracked: data/joined
Untracked: data/leafcutter/
Untracked: data/merged_combined_YL-SP-threeprimeseq.bg
Untracked: data/molPheno_noMP/
Untracked: data/mol_overlap/
Untracked: data/mol_pheno/
Untracked: data/nom_QTL/
Untracked: data/nom_QTL_opp/
Untracked: data/nom_QTL_trans/
Untracked: data/nuc6up/
Untracked: data/nuc_10up/
Untracked: data/other_qtls/
Untracked: data/pQTL_otherphen/
Untracked: data/peakPerRefSeqGene/
Untracked: data/perm_QTL/
Untracked: data/perm_QTL_GeneLocAnno_noMP_5percov/
Untracked: data/perm_QTL_GeneLocAnno_noMP_5percov_3UTR/
Untracked: data/perm_QTL_opp/
Untracked: data/perm_QTL_trans/
Untracked: data/perm_QTL_trans_filt/
Untracked: data/protAndAPAAndExplmRes.Rda
Untracked: data/protAndAPAlmRes.Rda
Untracked: data/protAndExpressionlmRes.Rda
Untracked: data/reads_mapped_three_prime_seq.csv
Untracked: data/smash.cov.results.bed
Untracked: data/smash.cov.results.csv
Untracked: data/smash.cov.results.txt
Untracked: data/smash_testregion/
Untracked: data/ssFC200.cov.bed
Untracked: data/temp.file1
Untracked: data/temp.file2
Untracked: data/temp.gencov.test.txt
Untracked: data/temp.gencov_zero.test.txt
Untracked: data/threePrimeSeqMetaData.csv
Untracked: data/threePrimeSeqMetaData55Ind.txt
Untracked: data/threePrimeSeqMetaData55Ind.xlsx
Untracked: data/threePrimeSeqMetaData55Ind_noDup.txt
Untracked: data/threePrimeSeqMetaData55Ind_noDup.xlsx
Untracked: data/threePrimeSeqMetaData55Ind_noDup_WASPMAP.txt
Untracked: data/threePrimeSeqMetaData55Ind_noDup_WASPMAP.xlsx
Untracked: output/picard/
Untracked: output/plots/
Untracked: output/qual.fig2.pdf
Unstaged changes:
Modified: analysis/28ind.peak.explore.Rmd
Modified: analysis/CompareLianoglouData.Rmd
Modified: analysis/NewPeakPostMP.Rmd
Modified: analysis/apaQTLoverlapGWAS.Rmd
Modified: analysis/cleanupdtseq.internalpriming.Rmd
Modified: analysis/coloc_apaQTLs_protQTLs.Rmd
Modified: analysis/dif.iso.usage.leafcutter.Rmd
Modified: analysis/diff_iso_pipeline.Rmd
Modified: analysis/explainpQTLs.Rmd
Modified: analysis/explore.filters.Rmd
Modified: analysis/flash2mash.Rmd
Modified: analysis/mispriming_approach.Rmd
Modified: analysis/overlapMolQTL.Rmd
Modified: analysis/overlapMolQTL.opposite.Rmd
Modified: analysis/overlap_qtls.Rmd
Modified: analysis/peakOverlap_oppstrand.Rmd
Modified: analysis/peakQCPPlots.Rmd
Modified: analysis/pheno.leaf.comb.Rmd
Modified: analysis/pipeline_55Ind.Rmd
Modified: analysis/swarmPlots_QTLs.Rmd
Modified: analysis/test.max2.Rmd
Modified: analysis/test.smash.Rmd
Modified: analysis/understandPeaks.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.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote
), click on the hyperlinks in the table below to view them.
File | Version | Author | Date | Message |
---|---|---|---|---|
html | b8d2b9f | Briana Mittleman | 2018-10-03 | Build site. |
Rmd | 2817554 | Briana Mittleman | 2018-10-03 | change window |
html | f8cf9e1 | Briana Mittleman | 2018-09-28 | Build site. |
Rmd | 94a7144 | Briana Mittleman | 2018-09-28 | all QTLs processed |
html | 28fef56 | Briana Mittleman | 2018-09-28 | Build site. |
Rmd | bfb5092 | Briana Mittleman | 2018-09-28 | 4su30qtl code |
html | e8acce6 | Briana Mittleman | 2018-09-28 | Build site. |
Rmd | 3914988 | Briana Mittleman | 2018-09-28 | 4su60qtl code |
html | a6038b1 | Briana Mittleman | 2018-09-27 | Build site. |
Rmd | 48be5d3 | Briana Mittleman | 2018-09-27 | add 4su processing |
Rmd | 79b6c46 | Briana Mittleman | 2018-09-27 | start to 4su |
html | d0019be | Briana Mittleman | 2018-09-27 | Build site. |
Rmd | 946b744 | Briana Mittleman | 2018-09-27 | add code for prot and RNA |
html | 47cd4bf | Briana Mittleman | 2018-09-26 | Build site. |
Rmd | fa7f707 | Briana Mittleman | 2018-09-26 | ribo QTL code |
I will use this analysis file to recall the other molecular QTLs using the same VCF files I am using for the APAqtls. This is important because we want to overlap QTLs called with the same genotype information.
processed (WASP+normalized) 4sU-seq (30m)
processed (WASP+normalized) 4sU-seq (60m)
processed (WASP+normalized) RNA-seq (Pickrell)
processed (WASP+normalized) RNA-seq (GEUVADIS)
processed (WASP+normalized) ribo-seq
LiftOver from (Battle et al., 2015) protein
I am download the processed data from http://eqtl.uchicago.edu/jointLCL/ and putting it in /project2/gilad/briana/threeprimeseq/data/molecular_phenos.
The protein file is already in the format needed for fastQTL. I need to change the headers to include the NA before the individuals.I will need to use:
export PATH=/project/gilad/software/midway1/qtltools-1.0:$PATH
bgzip phenotypes.bed && tabix -p bed phenotypes.bed.gz
To index the file for the program.
I will create a python script that adds the NA to the individuals.
def main(inF, outF):
infile= open(inF, "r")
fout = open(outF,'w')
for i, line in enumerate(infile):
if i == 0:
linelist=line.split()
for i, item in enumerate(linelist):
if i > 3:
linelist[i]="NA" + item
fout.write(" ".join(linelist) + '\n' )
else:
fout.write(line)
fout.close()
if __name__ == "__main__":
import sys
inF = sys.argv[1]
outF= sys.argv[2]
main(inF, outF)
Next step is to get the PCs to use as covariates in the analysis.
https://qtltools.github.io/qtltools/
This package is in /project/gilad/software/midway1/ and was installed by Peter Carbaneto from the RCC. I can add this to my path with:
export PATH=/project/gilad/software/midway1/qtltools-1.0:$PATH
I am going to use the QTLtools pca function. I need to run this on midway1.
QTLtools pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.bed.gz --scale --center --out /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.bed.PC.txt
#keep top 5 PCs for analysis
head -n 6 fastqtl_qqnorm_ribo_phase2.fixed.bed.PC.txt.pca > fastqtl_qqnorm_ribo_phase2.fixed.bed.5PCs.txt.pca
I then make a samples file wit the head of the PCA file. Remove 19192,19193 from sample file I need to make 1 vcf file with all of the chroms to run this.
riboQTL.nom.sh
#!/bin/bash
#SBATCH --job-name=riboQTL.nom
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=riboQTL.nom.out
#SBATCH --error=riboQTL.nom.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.bed.5PCs_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.noChr.bed.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_ribo_phase2.fixed.nominal.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/samples.txt
done
problem chr in pheno file and not in vcf
sed 's/^chr//' fastqtl_qqnorm_ribo_phase2.fixed.bed > fastqtl_qqnorm_ribo_phase2.fixed.noChr.bed
try changing /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.bed.5PCs.txt.pca first part of header to id like in the FastQTL site. and use tr to make it tap deliminated from " "
riboQTL.perm.sh
#!/bin/bash
#SBATCH --job-name=riboQTL.perm
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=riboQTL.perm.out
#SBATCH --error=riboQTL.perm.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.bed.5PCs_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_ribo_phase2.fixed.noChr.bed.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/perm/fastqtl_qqnorm_ribo_phase2.fixed.perm.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/samples.txt
done
I can process these to by generalizing my script from the apaQTLsAllInd.Rmd
I need to remove the Chr and add the NA.
fastqtl_qqnorm_RNAseq_phase2.txt
python addNAphenohead.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.txt /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.txt
sed 's/^chr//' /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.txt > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt
bgzip /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt
tabix /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt.gz
#midway1
QTLtools pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt.gz --scale --center --out /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.PC.txt
#keep top 5 PCs for analysis
head -n 6 fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.PC.txt.pca > fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.5PC.txt.pca
less fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.5PC.txt.pca | tr " " "\t" > fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.5PC_tab.txt.pca
I want to create a script that checks the samples in the file and then makes a samples.txt file if they are in the VCF.
First create a file with the names of the individuals in the vcf files. /project2/gilad/briana/YRI_geno_hg19/vcf.samples.txt
create_molQTLSamplefile.py
#infile is the pheno PCA file (because it is not zipped)
#outfile is the samples for that pheno if they are in the vcf samples
def main(inF, outF):
vcf_samp=open("/project2/gilad/briana/YRI_geno_hg19/vcf.samples.txt", "r")
for line in vcf_samp:
vcf_list=line.split()
infile= open(inF, "r")
fout = open(outF,'w')
for i, line in enumerate(infile):
if i ==0:
mol_samples=line.split()[1:]
final_samps=[]
for i in mol_samples:
if i in vcf_list:
final_samps.append(i)
fout.write("\n".join(final_samps))
if __name__ == "__main__":
import sys
inF = sys.argv[1]
outF= sys.argv[2]
main(inF, outF)
Run in the code dir:
python create_molQTLSamplefile.py ../data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.5PC_tab.txt.pca ../data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.Samples.txt
rnaQTL.nom.sh
#!/bin/bash
#SBATCH --job-name=rnaQTL.nom
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=rnaQTL.nom.out
#SBATCH --error=rnaQTL.nom.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.Samples.txt
done
rnaQTL.perm.sh
#!/bin/bash
#SBATCH --job-name=rnaQTL.perm
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=rnaQTL.perm.out
#SBATCH --error=rnaQTL.perm.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/perm/fastqtl_qqnorm_RNAseq_phase2.fixed.perm.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseq_phase2.Samples.txt
done
http://eqtl.uchicago.edu/jointLCL/fastqtl_qqnorm_prot.txt.gz
python /project2/gilad/briana/threeprimeseq/code/addNAphenohead.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.txt /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.txt
sed 's/^chr//' /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.txt > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.txt
bgzip /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.txt
tabix -p bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.txt.gz
#midway1
QTLtools pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.txt.gz --scale --center --out /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.PC.txt
#keep top 5 PCs for analysis
head -n 6 /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.PC.txt.pca > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.5PC.txt.pca
less /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.5PC.txt.pca | tr " " "\t" > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.5PC_tab.txt.pca
python create_molQTLSamplefile.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.5PC_tab.txt.pca /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.SAMP.txt
protQTL.nom.sh
#!/bin/bash
#SBATCH --job-name=protQTL.nom
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=protQTL.nom.out
#SBATCH --error=protQTL.nom.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_prot.fixed.nominal.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.SAMP.txt
done
protQTL.perm.sh
#!/bin/bash
#SBATCH --job-name=protQTL.perm
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=protQTL.perm.out
#SBATCH --error=protQTL.perm.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/perm/fastqtl_qqnorm_prot.fixed.nominal.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_prot.fixed.noChr.SAMP.txt
done
fastqtl_qqnorm_RNAseqGeuvadis_phase2.txt.gz
python addNAphenohead.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.txt /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.txt
sed 's/^chr//' /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.txt > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.txt
bgzip /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.txt
tabix -p bed -f /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.txt.gz
#midway1
QTLtools pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.txt.gz --scale --center --out /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.PC.txt
#keep top 5 PCs for analysis
head -n 6 /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.PC.txt.pca > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.5PC.txt.pca
less /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.5PC.txt.pca | tr " " "\t" > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.5PC_tab.txt.pca
python create_molQTLSamplefile.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.5PC_tab.txt.pca /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.SAMP.txt
RNAgQTL.nuc.sh
#!/bin/bash
#SBATCH --job-name=RNAgQTL.nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=RNAgQTL.nuc.out
#SBATCH --error=RNAgQTL.nuc.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_RNAseqGeuvadis.fixed.nominal.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.SAMP.txt
done
RNAgQTL.perm.sh
#!/bin/bash
#SBATCH --job-name=RNAgQTL.perm
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=RNAgQTL.perm.out
#SBATCH --error=RNAgQTL.perm.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/perm/fastqtl_qqnorm_RNAseqGeuvadis.fixed.perm.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_RNAseqGeuvadis_phase2.fixed.noChr.SAMP.txt
done
python addNAphenohead.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.txt /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.txt
sed 's/^chr//' /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.txt > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.txt
bgzip /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.txt
tabix -p bed -f /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.txt.gz
#midway1
QTLtools pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.txt.gz --scale --center --out /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.PC.txt
#keep top 5 PCs for analysis
head -n 6 /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.PC.txt.pca > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.5PC.txt.pca
less /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.5PC.txt.pca | tr " " "\t" > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.5PC_tab.txt.pca
python create_molQTLSamplefile.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.5PC_tab.txt.pca /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.SAMP.txt
4su60QTL_nom.sh
#!/bin/bash
#SBATCH --job-name=4su60QTL_nom
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=4su60QTL_nom.out
#SBATCH --error=4su60QTL_nom.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_4su60.fixed.nominal.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.SAMP.txt
done
4su60QTL_perm.sh
#!/bin/bash
#SBATCH --job-name=4su60QTL_perm
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=4su60QTL_perm.out
#SBATCH --error=4su60QTL_perm.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/perm/fastqtl_qqnorm_4su60.fixed.perm.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_60_phase2.fixed.noChr.SAMP.txt
done
python addNAphenohead.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.txt /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.txt
sed 's/^chr//' /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.txt > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.txt
bgzip /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.txt
tabix -p bed -f /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.txt.gz
#midway1
QTLtools pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.txt.gz --scale --center --out /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.PC.txt
#keep top 5 PCs for analysis
head -n 6 /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.PC.txt.pca > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.5PC.txt.pca
less /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.5PC.txt.pca | tr " " "\t" > /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.5PC_tab.txt.pca
python create_molQTLSamplefile.py /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.5PC_tab.txt.pca /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.SAMP.txt
4su30QTL_nom.sh
#!/bin/bash
#SBATCH --job-name=4su30QTL_nom
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=4su30QTL_nom.out
#SBATCH --error=4su30QTL_nom.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_4su30.fixed.nominal.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.SAMP.txt
done
4su30QTL_perm.sh
#!/bin/bash
#SBATCH --job-name=4su30QTL_perm
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=4su30QTL_perm.out
#SBATCH --error=4su30QTL_perm.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END
for i in $(seq 1 30)
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.5PC_tab.txt.pca --bed /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.txt.gz --out /project2/gilad/briana/threeprimeseq/data/molecular_QTLs/perm/fastqtl_qqnorm_4su30.fixed.perm.chunk$i.out --chunk $i 30 --window 5e5 --include-samples /project2/gilad/briana/threeprimeseq/data/molecular_phenos/fastqtl_qqnorm_4su_30_phase2.fixed.noChr.SAMP.txt
done
-cat the results and download them
sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS 10.14.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/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.2.0 Rcpp_0.12.19 digest_0.6.17 rprojroot_1.3-2
[5] backports_1.1.2 git2r_0.24.0 magrittr_1.5 evaluate_0.13
[9] stringi_1.2.4 fs_1.2.6 whisker_0.3-2 rmarkdown_1.11
[13] tools_3.5.1 stringr_1.4.0 glue_1.3.0 yaml_2.2.0
[17] compiler_3.5.1 htmltools_0.3.6 knitr_1.20