modified: MyMultiome/MultiomeATAC_mito.sh

modified: mitoConsensus/Finalize.sh modified: mitoConsensus/Preprocess.py
chenweng1991 · Oct 27, 2022 · 7ae909d · 7ae909d
1 parent ec30e47
commit 7ae909d
Showing 3 changed files with 69 additions and 16 deletions.
diff --git a/MyMultiome/MultiomeATAC_mito.sh b/MyMultiome/MultiomeATAC_mito.sh
@@ -1,15 +1,59 @@
 #!/bin/bash
 ## Note:  this version is for miseq or novaseq,  nextseq is slightly different
 
-echo "7 inputs: $1:Prefix, $2:Read1.fq, $3:Read2.fq, $4:i5.fq, $5:0; $6:threads, $7:MyMultiome path"
-name=$1
-Read1=$2
-Read2=$3
-ReadBarcode=$4
-Cut=$5 # Minimum uniq fragment per cell to be considered
-CORE=$6
-MyMultiome=$7
-bowtie2Index=$8
+Help()
+{
+  # Display Help
+  echo "This script trim and map the mtDNA fastqs, generating uniqmapped bam and QC plots "
+  echo
+  echo "MultiomeATAC_mito.sh -h for this page"
+  echo "Syntax: MultiomeATAC_mito.sh -n -1 -2 -i -c -t -m -b"
+  echo "Options"
+  echo "-n name: The prefix of all analyzed files"
+  echo "-1 Read1: Read1 of fastq file (150nt is recommended)"
+  echo "-2 Read2: Read2 of fastq file (150nt is recommended)"
+  echo "-i ReadBarcode: i5 of fastq file (24nt)"
+  echo "-c Cut: the cutoff the uniq fragment per cell"
+  echo "-t CORE: The number of cores to use"
+  echo "-m MyMultiome:The path to the folder of MyMultiome"
+  echo "-b bowtie2Index: the bowtie2 index path/prefix"
+  echo "-q quick, defult is false, if true then exit after uniqmapped.mito.bam, skipping QC step"
+}
+
+quick=0
+while getopts "hn:1:2:i:c:t:m:b:q" option; do
+  case $option in
+    h) # display help
+        Help
+        exit;;
+    n) # The prefix of all analyzed files
+        name=$OPTARG;;
+    1) # Read1 of fastq file
+        Read1=$OPTARG;;
+    2) # Read2 of fastq file
+        Read2=$OPTARG;;
+    i) # index5 fastq file
+        ReadBarcode=$OPTARG;;
+    c) # The cutoff the uniq fragment per cell
+        Cut=$OPTARG;;
+    t) # The number of cores to use
+        CORE=$OPTARG;;
+    m) # The path to the folder of MyMultiome
+        MyMultiome=$OPTARG;;
+    b) # The bowtie2 index path/prefix
+        bowtie2Index=$OPTARG;;
+    q) # if use this option then exit after uniqmapped.mito.bam, skipping QC step
+        quick=1;;
+   \?) # Invalid option
+        echo "Error: Invalid option"
+        exit;;
+  esac
+done
+
+## Exit if any of the necessary input is empty
+set -u
+: "$name$Read1$Read2$ReadBarcode$Cut$CORE$MyMultiome$bowtie2Index$quick" 
+
 
 ##Step 1 trim adaptor (Important)
 cutadapt --cores=$CORE -a CTGTCTCTTATA -A CTGTCTCTTATA -o $Read1.trim -p $Read2.trim $Read1 $Read2
@@ -37,6 +81,14 @@ echo "##Step5 Get uniq mapped bam and make bulk bigwig and call peaks"
 samtools index -@ $CORE $name.uniqmapped.bam
 samtools view -@ $CORE -b $name.uniqmapped.bam chrM > $name.uniqmapped.mito.bam
 
+if [[ quick -eq 1 ]]
+  then
+    echo "Skipping QC steps, exit"
+    exit
+  else
+    echo "Starting QC steps"
+fi
+
 ##Step6 Get raw bed file and Add cell barcode to the fragment bed files
 echo "##Step6 Get raw bed file and Add cell barcode to the fragment bed files"
 samtools sort -@ $CORE -n $name.uniqmapped.bam| bedtools bamtobed -bedpe -i stdin | awk -v OFS='\t' '{print $1,$2,$6,$7}' >$name.uniqmapped.RawBed

diff --git a/mitoConsensus/Finalize.sh b/mitoConsensus/Finalize.sh
@@ -1,7 +1,7 @@
 #!/bin/bash
 WD=$1
 Threads=$2
-CW_mgatk=/lab/solexa_weissman/cweng/Packages/CW_mgatk/
+mitoConsensus=$3
 
 if [ -z "$2" ]
 then
@@ -16,9 +16,9 @@ cat $WD/temp/sparse_matrices2.0/*RawGenotypes.VerySensitive >$WD/final/RawGenoty
 cat $WD/temp/sparse_matrices2.0/*RawGenotypes.Sensitive >$WD/final/RawGenotypes.Sensitive
 cat $WD/temp/sparse_matrices2.0/*RawGenotypes.Specific >$WD/final/RawGenotypes.Specific
 
-python3 $CW_mgatk/AddQualifiedTotalCts.py $WD
-Rscript $CW_mgatk/StrandBiases.R $WD 200 0.01
-python $CW_mgatk/RemoveStrandBias.py $WD
+python3 $mitoConsensus/AddQualifiedTotalCts.py $WD
+Rscript $mitoConsensus/StrandBiases.R $WD 200 0.01
+python $mitoConsensus/RemoveStrandBias.py $WD
 
 ##Final counting (Now it is included in the Finalize.sh script)
 let N=0
@@ -28,3 +28,4 @@ for i in $(seq 1 $Threads)
         let N=N+`samtools view $WD/temp/barcoded_bams/barcodes.$i.bam | wc -l`
         done
 echo -e $N > $WD/final/TotalRawBamRows.Mito
+
diff --git a/mitoConsensus/Preprocess.py b/mitoConsensus/Preprocess.py
@@ -8,7 +8,7 @@
 # barcodes="/lab/solexa_weissman/cweng/Projects/MitoTracing_Velocity/Run210706L1/Qualified_CellBC_ATAC_2k0.2.tb"
 output="TestFolder"
 mito_genome="rCRS"
-script_dir="/lab/solexa_weissman/cweng/Packages/CW_mgatk"
+#script_dir="/lab/solexa_weissman/cweng/Packages/CW_mgatk"
 
 @click.command()
 @click.version_option()
@@ -18,11 +18,11 @@
 @click.option('--output', '-o', default="mgatk_out", help='Output directory for analysis required for `call` and `bcall`. Default = mgatk_out')
 @click.option('--mito-genome', '-g', default = "rCRS", required=True, help='mitochondrial genome configuration. Choose hg19, hg38, mm10, (etc.) or a custom .fasta file; see documentation. Default = rCRS.')
 @click.option('--barcode-tag', '-bt', default = "BC", required=True, help='tag name see documentation. Default = BC.')
-
+@click.option('--script-dir', '-sd', default = "/lab/solexa_weissman/cweng/Packages/REDEEM-V/mitoConsensus/", required=True, help='')
 # For testing purpose
 
 
-def main(input,ncores,barcodes,output,mito_genome,barcode_tag):
+def main(input,ncores,barcodes,output,mito_genome,barcode_tag,script_dir):
     umi_barcode="XX"
     of = output; tf = of + "/temp"; bcbd = tf + "/barcoded_bams" # bcdb = barcoded bam directory
     folders = [of, tf, bcbd, of + "/final"] #,tf+"/sparse_matrices"