---
title: "R Notebook the process of hotspots detection in regulatory regions of the WCDT data "
output: html_notebook
author: Ruhollah
---

### Raw data
* __Annotated regions downloaded:__
    + 5kup, 2kdn, 5'UTR, 3'UTR from __ucsc_GenePred_GENCODEV28__ read/curated by `master_annotation_1.r` and written into `"master_annotation_list_reg.dat"`
    + CpGislands from __ucsc_Regulation_CpGIsland__ read/curated by `master_annotation_2_cpg.r` and appended into `"master_annotation_list_reg.dat"`
    + promoters, enhancers from __GeneHancer__ read/curated by `master_annotation_3_genehancer.r` and appended into `"master_annotation_list_reg.dat"`
    + genes from __Ensembl__ read/curated by `master_annotation_4_knowngenes.r` and appended into `"master_annotation_list_reg.dat"`
    + Data are then concatenated with the same number of features in columns
```{r}
# N.B. PLEASE DO NOT RUN WITHOUT TAKING BACKUPS OF THE RELEVENT FILES IN THE MAIN DIR. LOOK INTO THE SOURCE CODES FOR RELEVENT FILES

dir0 = 'C:/Users/ruhol/Box Sync/D/strelka_working/'
dir1 = '18.07.13.R--/'
source(paste(dir0,dir1,'0_sequential_run_masters_0.r', sep = ''))

```
  
  
  
  * __bash curation:__  
    +  `$ sort -V sorted_master_annotation > sorted_master_annotation.dat`
    *  `$ tail -n+2 sorted_master_annotation.dat  | sort -V | awk -F '[\t]' '{print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $6}' | uniq -c > uniqc_5col_sorted_mast_anot.dat`
    + `uniqc_5col_sorted_mast_anot.dat`  = final annotation file W/ ***CORRECTED*** 1-based for the UCSC data.
  
  * __Read & curate bash outputs:__
    + Add regreg_ID, remove very small & very large regregs,... to raw annotated columns of `uniqc_5col_sorted_mast_anot.dat` into an R datafram by `make_input_file_LASSO.R`'s folllowing excerpt :
  
```{r}
     
      raw_master_annot <- read.table('uniqc_5col_sorted_mast_anot.dat', sep = "")  
      colnames(raw_master_annot) <-  c("CHR", "start_reg", "end_reg", "strand", "regreg_type" )
      
      #raw_master_annot['regreg_type'] <- sub(":.*", "", raw_master_annot$regreg_type) # $sed -i 's/:.*//g' tttest.txt 

      #raw_master_annot['regreg_type'] <- str_replace_all(raw_master_annot$regreg_type ,"/","_") # $sed -i 's/\//_/g' tttest.txt 

      
      regreg_ID      <- paste(raw_master_annot$CHR, raw_master_annot$start_reg, raw_master_annot$end_reg, sep = "_")
      regreg_ID_type <- paste(raw_master_annot$CHR, raw_master_annot$start_reg, raw_master_annot$end_reg
                              , raw_master_annot$regreg_type,raw_master_annot$strand, sep = "_")
      raw_master_annot <- add_column(raw_master_annot, regreg_ID, .after = "regreg_type")
      raw_master_annot <- add_column(raw_master_annot, regreg_ID_type, .after = "regreg_ID")
  
      rownames(raw_master_annot) <- raw_master_annot$regreg_ID_type ###ALSO=THIS IS A TEST TO CATCH DUPLICATIONS IN 
        # INIT FILES AND RETURN TO CHECK THOSE FILES: this might give error for duplicated rownames 
        # => remove duplicated row from raw_master_annot
      
    
      
      ### write cleaned/uniq master annotation as the input file for the bedintersect w/ VCF files:
      # write.table(raw_master_annot[,c('CHR','start_reg', 'end_reg', 'strand', 'regreg_type')]
      #             , row.names = F, col.names = F, quote = F
      #             , file = "uniqc_5col_sorted_mast_anot.dat", sep = '\t') 
      #length regreg (1-based):
      raw_master_annot['reg_length'] <- raw_master_annot['end_reg'] - raw_master_annot['start_reg'] + 1
      
      rownames(raw_master_annot)<-c()
      
      ### truncate data based on the length of the regreg:
      
      trunc_mast_anot <- raw_master_annot[!raw_master_annot$reg_length<50,]
      trunc_mast_anot <- trunc_mast_anot[!trunc_mast_anot$reg_length>50000,]
      trunc2_mast_ant <- trunc_mast_anot[trunc_mast_anot$reg_length<10000,]
      
     
      
```
  
  * __One hot encoding__
    + All of the one hot encoding done on Avicenna server: `/avicenna/ruhollah/overlap`
    + Initially, the file `coverage_d_trunc2.dat` was created  by using the `bedtools coverage` over the entire genome. 
    + Then the information for each chr was extracted from `coverage_d_trunc2.dat` by `ruh_bash4.sh` and written to `chri.dat` `(i = 1:22,X,Y)`
    
     + Then files `chri.dat_overlap.dat` `(i = 1:22,X,Y)` for each chromosome created, showng the depth and the breadth of the overlaps.
     
     + Then files `chri_concated.dat` were created by the code `form_bed_concat.R`, showing the start and end of each partial regulatory region on `chri`. These files show the start, end, length, breadth, scand depth (= total number of distincts regions controbuting in this partial region). 
     + Then files `chri.dat_overlap.txt`  were created from `chri_concated.dat`  In these files, the end of partial regreg are corrected; and also the partial regions w/ breadth=1 were removed. 
    + Then  tab-delimited files `chri_chopped.dat1` were created out of the latter files via  `make_chopped_trunc2_mast_anot.sh`.
    + This latter bash script also concatenate and `sort -V` all `.dat1` into `chopped_trunc2_mast_anot.dat`. 
    Finally, the recent bash script does the `bedintersect` b/w  `chopped_trunc2_mast_anot.dat` & the previously-created annotation file to assign the name of the each chopped (=partial) regulatory region to it in the last column:
    `$bedtools intersect -a chopped_trunc2_mast_anot.dat -b trunc2_uniq_mast_ant.dat -wo > bedintersect_chopped_trunc2_mast_ant.dat`

    + Finally, from this recent `.dat` file, the final file `merge4col_bedint_chopped_trunc2.dat` was created in bash by removing extra columns and merging the repeated lines. An extra column was added to the latter file to give the final file of `chopped_trunc2_mast_anot_pregnames.dat`. This column shows the number of distinct original regregs overlapping in each partial, chopped regreg. 
   
   
     
  * __Epigenetic features: download & curation:__
    + All feature types from GRCh38 PC3 cell line from the ENCODE project were downloaded as `bed` files, __supplemented__ by the corresponding dataset names and target feature names via the bash script `strelka_working\Encode_data_peaks\PC3_bed.narrowPeaks\GRCh38Only\ruh_bash2.sh`. 
    
    + Then these files were __sorted,  briefed (only 5 columns), & cleaned (got rid of unknown chromosomes)__. Then these files were __concatenated and sorted__ to give--> `sorted_ALL_brief_clean.bed`. Then this file was __briefed__ by removing the column of the filename to give--> `4cols_sorted_ALL_brief_clean.bed`. Then __duplicated lines were removed__ from this file to give --> `norep_sorted_All_breif_clean.bed`.
  
  
      <!-- 
      + This latter file then was __`bedintersected`__ by the `VCF` files to give `bedintersect_narrowpeaks.out`. 
      + This file was then __curated__ in `R` to remove repeated bedintersected lines and was writtern as `5cols_bedintersect_norepMasterAnnot_norepNarrowpeaks_PC3.out`, which was then `$ uniq -c > uniqc_5cols_bedintersect_norepMasterAnnot_norepNarrowpeaks_PC3.out` 
      -->

  
  * __bedintersect:__ 
  
  
    <!-- This is for the original:
    `bedintersect`, before one hot encoding:  `uniqc_5col_sorted_mast_anot.dat` bedintersected with *all* of  `VCF` files *at once* by the bashscript `ruh_bash1_strelka_bedtools_bedintersect.sh`.  
    -->

   + `chopped_trunc2_mast_anot_pregnames.dat` bedintersected with *all* of  `VCF` files *at once* by the bashscript `ruh_bash1_strelka_bedtools_bedintersect.sh` to give--> `bedintersect_chopped_trunc2_w_VCF_files.out`
   
   + Similarly, epigenetic features were bedinteresected with chopped regulatory regions. 

  
  * __Read & curate bedintersect results in `BedIntresect_analysis.R`:__ 
    +  Filter out  rows w/ AF1000G >= 0.01  & SNPs
     + Saved as `res2_bedintersect_chopped_regreg_trunc2.RData`(this code does other things as well, like subtype classification of mCRPC and computing the related correlation maps)
     + Output was written into  `res2_4col_bedintersect_chopped_trunc2.dat` >  `bash-sort-V sortV_res2_4col_bedintersect_chopped_trunc2.dat` > `bash-uniq-c uniqc_sortV_res2_4col_bedintersect_chopped_trunc2.dat` to give the numebr of mutation (`nmt`) per partial(=chopped) regreg 
     
  
```{r}
### N.B. MAKE SURE TO CHNAGE THE CODE AT BEGINNING OF THE BedIntresect_analysis.R TO CHOOSE THE RIGHT INPUT FILE, AS THIS CODE IS GENERAL AND IS USED TO ANALYZE DIFFERENT BEDINTERSECT RESULTs 

source(paste(dir0,'BedIntresect_analysis.R' , sep = ''))
```
 
  * __GC content per chopped (=partial) regreg (i.e., one hot encoded):__
	
    + The fasta sequence of each chopped (=partial) regreg was computed based on hg38 genome by the bash twoBitToFa at avicenna as: 
    + `$ /usr/local/miniconda2/bin/twoBitToFa -noMask -bed=chopped_4col_trunc2_mast_ant.dat /avicenna/genomes/hg38/hg38.2bit seq_chopped_trunc2_mst_ant.out`

    + GC% calculated by `/avicenna/ruhollah/overlap/dat1files/gc_content.pl` on the fasta sequences of each chopped (=partial) regreg as `$ ./gc_content.pl seq_chopped_trunc2_mst_ant.out` 
       
  * __Dinucleotide distribution:__    
    + `kmer-counter` & `jellyfish-2.2.10` --> dinucleotide density 

  * __Some further curation:__
    + Remove segments (=partial regreg) w/ only one-type bp (AAAA.. or TTT…) (a total of 42 segements)
    + Remove all chopped segments below 10bp 
    + Data were saved as `x_y_preg.RData`
  
  

### Model building
#### Data exploratory:
  * __Correlation matrices:__
    + For epigenetic features, dinucleotide & GC%, and one-hot encoded regregs. 
    + Done with the following exceprt of the `make_input_file_LASSO.R`:
```{r}
  cormat.epi <- rcorr(as.matrix(xy_preg[,c(22:38)])) ##findng pariwise correlation b/w all epi-features |--> symmetric matrix of corr.coeff
  arrind <- which(cormat.epi[[1]] >.7 & lower.tri(cormat.epi[[1]] ), arr.ind = T) ##extracting highly correlated epi-features (r>.5, .7 , ...), excluding the self-correlated
  a <- cbind(row.names(arrind), colnames(cormat.epi[[1]])[arrind[,2]])  # forming the matrix of the names of highly correlated epi-features. 
  
  
  cormat.dinuc <- rcorr(as.matrix(xy_preg[,c(10,39:54)])) ##GCp , AA-TT
  arrind <- which(cormat.dinuc[[1]] >.5 & lower.tri(cormat.dinuc[[1]] ), arr.ind = T)
  a <- cbind(row.names(arrind), colnames(cormat.dinuc[[1]])[arrind[,2]]) 
  
  cormat.regreg <- rcorr(as.matrix(xy_preg[,c(11:19)]))
  arrind <- which(abs(cormat.regreg[[1]]) >.35 & lower.tri(cormat.regreg[[1]] ), arr.ind = T)
  a <- cbind(row.names(arrind), colnames(cormat.regreg[[1]])[arrind[,2]]) 
```
  
  * __PCA of dinucleotides, GC%, and Epigenetic features__
    + Conducted in `make_input_file_LASSO.R`. 
    + Out of 17 dinucleotide & GC%, the first 7 PCs were selected that capture >80% variations
    + Out of 16 epigenetic features, the first 10 PCs """""""""""""""""""""""""""""""""""""""
    + Results were saved as `xy_pca.RData` 
    
#### Model fitting

  * __Fitting the model:__ 
    + Finding the best distribution for the response, `dnmt`  (density of mutations in partial regreg), by trying transformation using different concave functions and evaluate the results through Cullen&Frey graph
    + `dnmt^0.1` gives the best fit for `gamma` distribution 
    + Ignore the `dnmt = 0` values (see `make_input_file_LASSO.R`): 
```{r}
glm_pc_dnmt.1_Gamma_no0 <- glm(dnmt.1 ~ Gene+Enhancer+tkdns+fpUTR+fkups+Promoter+CpG+thpUTR+PromHancer+nn_PC1_dinuc+nn_PC2_dinuc+nn_PC3_dinuc+nn_PC4_dinuc+nn_PC5_dinuc+nn_PC6_dinuc+nn_PC7_dinuc+nn_PC1_epi+nn_PC2_epi+nn_PC3_epi+nn_PC4_epi+nn_PC5_epi+nn_PC6_epi+nn_PC7_epi+nn_PC8_epi+nn_PC9_epi+nn_PC10_epi
                           , family = Gamma(), data = xy_pc[xy_pc$nmt>0,])
```
     
      Results saved as `glm_brief.RData`
      Use `ruh_outlier` function to sort outliers from `glm` output:
    
```{r}
finaloutliers <- ruh_outlier(origxy = xy_pc, expdnmt = .1, glm_result = glm_pc_dnmt.1_Gamma_no0, minnmt = 3,niq = 2.5)
```
  
    + `expdnmt` = the exponent of the concave transformation function used to fit `dnmt` into appropriate distribution
    + `minnmt` = minimum  number of mutations in the partial regreg that must be imposed
    + `niq` = Determines the  interquatile range of the residuals 
    
    
### Gene expression analysis
```{r}
source(paste(dir0, 'assigning_genes_and_expressions_to_outliers-NOtebook.Rmd', sep = '')
```

   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
  