updates

Author: p4rkerw

allele_specific_analysis/step1a_gatk_rna_genotype.sh

@@ -34,7 +34,7 @@
 # -v g/reference/gatk:$HOME/gatk_bundle \
 # -v $SCRATCH1:$SCRATCH1 \
 # -e SCRATCH1="/g/scratch" \
-# -it p4rkerw/gatk:4.1.8.1
+# --rm -it p4rkerw/gatk:4.1.8.1
 
 # make sure that miniconda is in the path
 export PATH=$PATH:/gatk:/opt/miniconda/envs/gatk/bin:/gatk:/opt/miniconda/envs/gatk/bin:/opt/miniconda/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin
@@ -62,7 +62,7 @@ gatk MarkDuplicates \
       --CREATE_INDEX true
 
 # limit to specified intervals
-# this will speed up analysis be eliminating alt contigs
+# this will speed up analysis by eliminating alt contigs
 CONTIGS=(chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY)      
 printf '%s\n' "${CONTIGS[@]}" > /tmp/interval.list
 
@@ -93,18 +93,37 @@ gatk ApplyBQSR \
 
 # single sample gvcf variant calling 
 # annotate the sites with dbsnp
-gatk HaplotypeCaller  \
-   -R ref/genome.fa \
-   -I $outdir/bqsr.bam \
-   -O $outdir/output.g.vcf.gz \
-   -ERC GVCF \
-   -D gatk_bundle/resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf
+PROCESSES=4  # Number of threads you want
+for INTERVAL in ${CONTIGS[*]}; do
+echo $INTERVAL
+   gatk HaplotypeCaller  \
+      -R ref/genome.fa \
+      -I $outdir/bqsr.bam \
+      -O $outdir/output.g.$INTERVAL.vcf.gz \
+      -ERC GVCF \
+      -D gatk_bundle/resources_broad_hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf \
+      -L $INTERVAL \
+      --verbosity ERROR &
+      while (( $(jobs |wc -l) >= (( ${PROCESSES} + 1 )) )); do
+         sleep 0.1
+      done
+done
+
+# merge vcfs from parallel haplotypecaller 
+(ls $outdir/output.g.chr*.vcf.gz) > /tmp/vcf.list
+gatk GatherVcfs \
+   -I /tmp/vcf.list \
+   -O $outdir/output.g.vcf.gz 
+
+# index merged vcf
+gatk IndexFeatureFile \
+   -I $outdir/output.g.vcf.gz 
 
 # single sample genotyping
 gatk GenotypeGVCFs \
    -R ref/genome.fa \
    -V $outdir/output.g.vcf.gz \
-   -O $outdir/output.vcf.gz
+   -O $outdir/output.vcf.gz 
 
 # perform hard filtering using the qual-by-depth QD score and window for snp clustering
 gatk VariantFiltration \
@@ -116,6 +135,6 @@ gatk VariantFiltration \
    --filter "FS > 30.0" \
    --filter-name "QD" \
    --filter "QD < 2.0" \
-   -O $outdir/filter.output.vcf.gz
+   -O $outdir/$SAMPLE_NAME.filter.output.vcf.gz
 
 

analysis/corr_chromVar_TF_exp.R

@@ -134,6 +134,16 @@ df.final <- dplyr::distinct(df.final, celltype, gene_motif, .keep_all = TRUE)
 df.sig <- dplyr::filter(df.final, pval < 0.05) %>%
   arrange(-num_celltypes, -corr)
 
+df.final$celltype <- as.factor(df.final$celltype)
+levels(df.final$celltype) <- c("PT","PTVCAM1","PEC","TAL","DCT1","DCT2","CNT","PC","ICA","ICB",
+                               "MES","FIB","ENDO","PODO","LEUK")
+
+
+df.sig$celltype <- as.factor(df.sig$celltype)
+levels(df.sig$celltype) <- c("PT","PTVCAM1","PEC","TAL","DCT1","DCT2","CNT","PC","ICA","ICB",
+                               "MES","FIB","ENDO","PODO","LEUK")
+
+
 # calculate pearson r2 for all tf-gene combos
 pearson <- cor.test(df.sig$chromvar, df.sig$rna, method="pearson", conf.level=0.95)
 max_chrom <- max(abs(df.sig$max_chrom)) + 1
@@ -164,7 +174,7 @@ pearson <- cor.test(df.pos$chromvar, df.pos$rna, method="pearson", conf.level=0.
 max_chrom <- max(abs(df.pos$max_chrom)) + 1
 max_exp <- max(abs(df.pos$max_exp)) + 1
 
-p1 <- ggplot(df.pos, aes(x=chromvar, y=rna)) +
+p2 <- ggplot(df.pos, aes(x=chromvar, y=rna)) +
   geom_smooth(method="lm", color = "darkgray") +
   geom_point(aes(color=celltype)) +
   xlab("chromVAR activity (avg_logFC)") +
@@ -176,7 +186,7 @@ p1 <- ggplot(df.pos, aes(x=chromvar, y=rna)) +
   theme_minimal() +
   geom_hline(yintercept=0) +
   geom_vline(xintercept=0) 
-p1
+p2
 
 
 
@@ -186,7 +196,7 @@ pearson <- cor.test(df.neg$chromvar, df.neg$rna, method="pearson", conf.level=0.
 max_chrom <- max(abs(df.neg$max_chrom)) + 1
 max_exp <- max(abs(df.neg$max_exp)) + 1
 
-p2 <- ggplot(df.neg, aes(x=chromvar, y=rna)) +
+p3 <- ggplot(df.neg, aes(x=chromvar, y=rna)) +
   geom_smooth(method="lm", color = "darkgray") +
   geom_point(aes(color=celltype)) +
   xlab("chromVAR activity (avg_logFC)") +
@@ -198,13 +208,31 @@ p2 <- ggplot(df.neg, aes(x=chromvar, y=rna)) +
   theme_minimal() +
   geom_hline(yintercept=0) +
   geom_vline(xintercept=0) 
-p2
+p3
 
+gene_motif <- df.final$gene_motif[grepl("NR3C1", df.final$gene_motif)]
+toplot <- dplyr::filter(df.final, motif == "MA0113.3")
+max_chrom <- unique(toplot$max_chrom)
+max_exp <- unique(toplot$max_exp)
+p4 <- ggplot(toplot, aes(x=chromvar, y=rna)) +
+  geom_smooth(method="lm", color = "darkgray") +
+  geom_point(aes(color=celltype)) +
+  geom_text_repel(aes(label=celltype)) +
+  xlab("chromVAR activity (avg_logFC)") +
+  ylab("Gene expression (avg_logFC)") +
+  ggtitle(unique(toplot$gene), subtitle=paste0("Pearson r^2=",unique(toplot$cor)," pval=",unique(toplot$pval))) +
+  xlim(c(-max_chrom,max_chrom)) +
+  ylim(c(-max_exp,max_exp)) +
+  theme_minimal() +
+  geom_hline(yintercept=0) +
+  geom_vline(xintercept=0) 
+p4
 
-toplot <- dplyr::filter(df.final, gene_motif == "ZEB1_MA0103.3")
+gene_motif <- df.final$gene_motif[grepl("NR3C2", df.final$gene_motif)]
+toplot <- dplyr::filter(df.final, motif == "MA0727.1")
 max_chrom <- unique(toplot$max_chrom)
 max_exp <- unique(toplot$max_exp)
-p1 <- ggplot(toplot, aes(x=chromvar, y=rna)) +
+p5 <- ggplot(toplot, aes(x=chromvar, y=rna)) +
   geom_smooth(method="lm", color = "darkgray") +
   geom_point(aes(color=celltype)) +
   geom_text_repel(aes(label=celltype)) +
@@ -216,5 +244,5 @@ p1 <- ggplot(toplot, aes(x=chromvar, y=rna)) +
   theme_minimal() +
   geom_hline(yintercept=0) +
   geom_vline(xintercept=0) 
-p1
+p5