From 68a092037bbcbe3dde940c42bf8d025edc41cc2d Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Mon, 17 Jan 2022 19:11:10 +0000
Subject: [PATCH] finding seq discrepancy in MSA for embb

---
 scripts/functions/logoP_snp.R             |  6 +++---
 scripts/functions/tests/test_logo_plots.R | 25 ++++++++++++++++++++++-
 scripts/gene_targets_names.txt            |  3 +--
 scripts/run_mutate2.sh                    | 11 +++++++---
 4 files changed, 36 insertions(+), 9 deletions(-)

diff --git a/scripts/functions/logoP_snp.R b/scripts/functions/logoP_snp.R
index fa01541..bee8e9f 100644
--- a/scripts/functions/logoP_snp.R
+++ b/scripts/functions/logoP_snp.R
@@ -195,10 +195,10 @@ LogoPlotSnps <- function(plot_df
   cat('\nDone: p0')
   
   # further customisation
-  mut_logo_p <<- p0 + theme(legend.position = leg_pos
+  mut_logo_p = p0 + theme(legend.position = leg_pos
                           , legend.direction = leg_dir
                           #, legend.title = element_blank()
-                          , legend.title = element_text(size = y_tts
+                          , legend.title = element_text(size = leg_tts
                                                         , colour = ytt_col)
                           , legend.text = element_text(size = leg_ts)
 
@@ -246,7 +246,7 @@ LogoPlotSnps <- function(plot_df
   cat('\nDone: p1')
   
   # further customisation
-  wt_logo_p <<- p1 + 
+  wt_logo_p = p1 + 
     
     theme(legend.position = "none"
           , legend.direction = leg_dir
diff --git a/scripts/functions/tests/test_logo_plots.R b/scripts/functions/tests/test_logo_plots.R
index 9065394..43ce22d 100644
--- a/scripts/functions/tests/test_logo_plots.R
+++ b/scripts/functions/tests/test_logo_plots.R
@@ -1,4 +1,6 @@
-source("~/git/LSHTM_analysis/config/gid.R")
+#source("~/git/LSHTM_analysis/config/gid.R")
+source("~/git/LSHTM_analysis/config/alr.R")
+
 source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
 
 ################################
@@ -56,3 +58,24 @@ LogoPlotSnps(plot_df = merged_df3
              , leg_tts = 16 # leg title size
 )
 
+
+
+########################################
+# Logo plot MSA
+# Mutant and wild-type
+# wild-type and mutant aa
+# script: logoP_msa.R
+########################################
+# msa1 = read.csv("/home/tanu/git/Data/cycloserine/output/alr_msa.csv", header = F)
+# head(msa1)
+# msa_seq= msa1$V1
+# head(msa_seq)
+# 
+# msa2 = read.csv("/home/tanu/git/Data/cycloserine/input/alr.1fasta", header = F)
+# head(msa2)
+# wt_seq = msa2$V1
+# head(wt_seq)
+# 
+# # BOTH WORK
+# LogoPlotMSA(msa_seq, wt_seq)
+# LogoPlotMSA(msa1, msa2)
diff --git a/scripts/gene_targets_names.txt b/scripts/gene_targets_names.txt
index 9b9a970..25c6e02 100644
--- a/scripts/gene_targets_names.txt
+++ b/scripts/gene_targets_names.txt
@@ -1,7 +1,6 @@
-embb ethambutol
 rpob rifampicin
 alr cycloserine
 katg isoniazid
 pnca pyrazinamide
 gid streptomycin
-
+embb ethambutol
diff --git a/scripts/run_mutate2.sh b/scripts/run_mutate2.sh
index 5dd5c60..2038a11 100755
--- a/scripts/run_mutate2.sh
+++ b/scripts/run_mutate2.sh
@@ -30,6 +30,8 @@ while read -r gene drug; do
   echo "Running mutate.py on data file $MSA_MAP"
   python3 mutate.py -v -o ${DATA_DIR}/${drug}/output/${gene}_msa_interim.csv $MSA_MAP $DATA_DIR/${drug}/input/${gene}_f2.fasta
   echo "mutate.py completed"
+  sed -E 's/>.*//g;/^$/d' ${DATA_DIR}/${drug}/output/${gene}_msa_interim.csv > ${DATA_DIR}/${drug}/output/${gene}_msa.csv
+  wc -l ${DATA_DIR}/${drug}/output/${gene}_msa.csv
   echo
 
 done < gene_targets_names.txt
@@ -37,14 +39,17 @@ done < gene_targets_names.txt
 # Stop here so we don't run the examples below :)
 exit
 
+########################################################################
+#
+########################################################################
 # make sure there is no new line at the end of the mutation file (snps.csv)
 # check
-cat  output/gid_metadata.csv | rev| cut -d, -f1 |rev | tail -n +2 |sort | head
+cat  output/gid_metadata.csv | rev | cut -d, -f1 |rev | tail -n +2 |sort | head
 
-cat  output/gid_metadata.csv | rev| cut -d, -f1 |rev | tail -n +2 |sort | uniq -c  > output/gid_metadata_mut_count.csv
+cat  output/gid_metadata.csv | rev | cut -d, -f1 |rev | tail -n +2 |sort | uniq -c  > output/gid_metadata_mut_count.csv
 
 
-cat  output/gid_metadata.csv | rev| cut -d, -f1 |rev | tail -n +2 |sort > gid_msa_snp.csv
+cat  output/gid_metadata.csv | rev | cut -d, -f1 |rev | tail -n +2 |sort > gid_msa_snp.csv
 sed -i 's/^/gid,/' gid_msa_snp.csv 
 
 #cp gid_msa_snp.csv gid_mut_map.csv