Last updated: 2018-10-03
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: 2817554
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: KalistoAbundance18486.txt
Untracked: analysis/genometrack_figs.Rmd
Untracked: analysis/ncbiRefSeq_sm.sort.mRNA.bed
Untracked: analysis/snake.config.notes.Rmd
Untracked: analysis/verifyBAM.Rmd
Untracked: data/18486.genecov.txt
Untracked: data/APApeaksYL.total.inbrain.bed
Untracked: data/NuclearApaQTLs.txt
Untracked: data/RNAkalisto/
Untracked: data/TotalApaQTLs.txt
Untracked: data/Totalpeaks_filtered_clean.bed
Untracked: data/YL-SP-18486-T-combined-genecov.txt
Untracked: data/YL-SP-18486-T_S9_R1_001-genecov.txt
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/ensemble_to_genename.txt
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/nom_QTL/
Untracked: data/nom_QTL_opp/
Untracked: data/nom_QTL_trans/
Untracked: data/nuc6up/
Untracked: data/other_qtls/
Untracked: data/peakPerRefSeqGene/
Untracked: data/perm_QTL/
Untracked: data/perm_QTL_opp/
Untracked: data/perm_QTL_trans/
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: output/picard/
Untracked: output/plots/
Untracked: output/qual.fig2.pdf
Unstaged changes:
Modified: analysis/28ind.peak.explore.Rmd
Modified: analysis/39indQC.Rmd
Modified: analysis/PeakToGeneAssignment.Rmd
Modified: analysis/cleanupdtseq.internalpriming.Rmd
Modified: analysis/dif.iso.usage.leafcutter.Rmd
Modified: analysis/diff_iso_pipeline.Rmd
Modified: analysis/explore.filters.Rmd
Modified: analysis/overlapMolQTL.Rmd
Modified: analysis/overlap_qtls.Rmd
Modified: analysis/peakOverlap_oppstrand.Rmd
Modified: analysis/pheno.leaf.comb.Rmd
Modified: analysis/test.max2.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 | 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 Sierra 10.12.6
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.1.1 Rcpp_0.12.18 digest_0.6.16
[4] rprojroot_1.3-2 R.methodsS3_1.7.1 backports_1.1.2
[7] git2r_0.23.0 magrittr_1.5 evaluate_0.11
[10] stringi_1.2.4 whisker_0.3-2 R.oo_1.22.0
[13] R.utils_2.7.0 rmarkdown_1.10 tools_3.5.1
[16] stringr_1.3.1 yaml_2.2.0 compiler_3.5.1
[19] htmltools_0.3.6 knitr_1.20
This reproducible R Markdown analysis was created with workflowr 1.1.1