05-1-Partition_CMRG_sequence_dipcall.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "97a08c22",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing /code/pgrtk-0.3.1-cp38-cp38-linux_x86_64.whl\n",
      "Installing collected packages: pgrtk\n",
      "  Attempting uninstall: pgrtk\n",
      "    Found existing installation: pgrtk 0.3.1\n",
      "    Uninstalling pgrtk-0.3.1:\n",
      "      Successfully uninstalled pgrtk-0.3.1\n",
      "Successfully installed pgrtk-0.3.1\n"
     ]
    }
   ],
   "source": [
    "!pip install --force /code/pgrtk-0.3.2-cp38-cp38-linux_x86_64.whl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0eee28ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pgrtk\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import SVG\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9c7b1220",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'pgrtk 0.3.1 (docs:e7ec376, release build, linux [x86_64] [rustc 1.62.0 (a8314ef7d 2022-06-27)])'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pgrtk.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e48093b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"/results\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "dc52796e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 47s, sys: 29.6 s, total: 4min 17s\n",
      "Wall time: 25.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "sdb = pgrtk.SeqIndexDB()\n",
    "sdb.load_from_fastx(\"/results/output_dir/out_seqs_with_padding.fa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cacb1880",
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_info = sdb.seq_info.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a9c0a813",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "37638"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(seq_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ecd75a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "4e6cc0a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_to_seqeunces = {}\n",
    "for sid in seq_info:\n",
    "    ctg, src, _ = seq_info[sid]\n",
    "    #print(ctg)\n",
    "    gene, seq_name = ctg.split(\"_\",1)\n",
    "    #if gene not in (\"SNTG2\", \"SLC27A5\"):\n",
    "    #    continue\n",
    "    try:\n",
    "        sample, haplotype, _ = seq_name.split(\"#\", 2)\n",
    "        haplotype = int(haplotype) - 1\n",
    "        sample_to_seqeunces.setdefault(sample, [[],[]])\n",
    "        seq = pgrtk.u8_to_string(sdb.get_seq(src, ctg))\n",
    "        sample_to_seqeunces[sample][haplotype].append( (ctg, seq) )\n",
    "        #print(sample, haplotype)\n",
    "    except:\n",
    "        pass\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "5e1c84b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir -p /results/output2/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "7135a81a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for sample in sample_to_seqeunces:\n",
    "    for haplotype in (0, 1):\n",
    "        with open(\"/results/output2/{}_{}_CMRG.fa\".format(sample, haplotype), \"w\") as f:\n",
    "            for ctg, seq in sample_to_seqeunces[sample][haplotype]:\n",
    "                print(f\">{ctg}\", file=f)\n",
    "                print(f\"{seq}\", file=f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5960b4d",
   "metadata": {},
   "source": [
    "Install `dipcall` with `conda install -c bioconda dipcall`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "941768bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "!/code/agc-1.1_x64-linux/agc getset /data/grch38.agc  GCA_000001405.15_GRCh38_no_alt_analysis_set.fna > /results/grch38.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bb8b72b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /results/output2/\n",
    "minimap2 -d grch38.mmidx ../grch38.fa \n",
    "samtools faidx ../grch38.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "d5f82458",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting /results/output2/run_dipcall\n"
     ]
    }
   ],
   "source": [
    "%%writefile /results/output2/run_dipcall\n",
    "#!/usr/bin/env perl\n",
    "\n",
    "use strict;\n",
    "use warnings;\n",
    "use Getopt::Std;\n",
    "\n",
    "my $version = \"0.3\";\n",
    "\n",
    "my %opts = (t=>8);\n",
    "getopts(\"t:d:x:muz:W:a\", \\%opts);\n",
    "\n",
    "die(\"Usage: run-dipcall [options] <prefix> <ref.fa> <pat.fa> <mat.fa>\n",
    "Options:\n",
    "  -t INT     number of threads [$opts{t}]\n",
    "  -d FILE    unimap/minimap2 index for ref.fa []\n",
    "  -a         call on all contigs regardless of naming\n",
    "  -x FILE    PAR on chrX; assuming male\n",
    "  -z INT     Z-drop [mapper default]\n",
    "  -m         use minimap2 for mapping (default)\n",
    "  -u         use unimap for mapping\n",
    "  -W FILE    repetitive k-mer list; use winnowmap for mapping\n",
    "\") if @ARGV < 4;\n",
    "\n",
    "# test file existence\n",
    "my $pre = shift(@ARGV);\n",
    "my $ref = shift(@ARGV);\n",
    "my @hap = @ARGV;\n",
    "die \"ERROR: failed to read the reference file '$ref'\\n\" unless -f $ref;\n",
    "die \"ERROR: please index the reference with 'samtools faidx'\\n\" unless -f \"$ref.fai\";\n",
    "die \"ERROR: failed to read the 1st haplotype '$hap[0]'\\n\" unless -f $hap[0];\n",
    "die \"ERROR: failed to read the 2nd haplotype '$hap[1]'\\n\" unless -f $hap[1];\n",
    "my $is_male = defined($opts{x})? 1 : 0;\n",
    "if (defined $opts{x}) {\n",
    "\tdie \"ERROR: failed to read PAR\\n\" unless -f $opts{x};\n",
    "\t$is_male = 1;\n",
    "}\n",
    "\n",
    "# find the root directory\n",
    "my $exepath = $0 =~/^\\S+\\/[^\\/\\s]+/? $0 : &which($0);\n",
    "#my $root = $0 =~/^(\\S+)\\/[^\\/\\s]+/? $1 : undef;\n",
    "#$root = $exepath =~/^(\\S+)\\/[^\\/\\s]+/? $1 : undef if !defined($root);\n",
    "my $root = \"/opt/conda/bin/\";\n",
    "die \"ERROR: failed to locate the root directory\\n\" if !defined($root);\n",
    "# mapper settings\n",
    "my $mapper = \"minimap2\";\n",
    "my $mm2_opt = q/-xasm5 --cs -t$(N_THREADS) -z200000,10000/;\n",
    "if (defined $opts{u}) {\n",
    "\t$mapper = \"unimap\";\n",
    "} elsif (defined $opts{W}) {\n",
    "\t$mapper = \"winnowmap\";\n",
    "\t$mm2_opt .= \" -r2k -W $opts{W}\";\n",
    "}\n",
    "$mapper = \"minimap2\" if defined($opts{m});\n",
    "$mm2_opt .= \" -z$opts{z},200\" if defined($opts{z});\n",
    "my $mm2_idx = defined($opts{d})? $opts{d} : $ref;\n",
    "\n",
    "# write Makefile\n",
    "my @mak = ();\n",
    "\n",
    "push(@mak, \"ROOT=$root\");\n",
    "push(@mak, \"N_THREADS=$opts{t}\");\n",
    "push(@mak, \"MAPPER=$mapper\");\n",
    "push(@mak, \"MM2_IDX=$mm2_idx\");\n",
    "push(@mak, \"MM2_OPT=$mm2_opt\");\n",
    "push(@mak, \"REF_FA=$ref\");\n",
    "push(@mak, \"\");\n",
    "\n",
    "push(@mak, \"all:$pre.dip.bed $pre.dip.vcf.gz\", \"\");\n",
    "\n",
    "push(@mak, \"$pre.hap1.paf.gz:$hap[0]\");\n",
    "push(@mak, q{\t$(ROOT)/$(MAPPER) -c --paf-no-hit $(MM2_OPT) $(MM2_IDX) $< 2> $@.log | gzip > $@});\n",
    "push(@mak, \"$pre.hap2.paf.gz:$hap[1]\");\n",
    "push(@mak, q{\t$(ROOT)/$(MAPPER) -c --paf-no-hit $(MM2_OPT) $(MM2_IDX) $< 2> $@.log | gzip > $@});\n",
    "push(@mak, \"\");\n",
    "\n",
    "push(@mak, \"$pre.hap1.sam.gz:$hap[0]\");\n",
    "push(@mak, q{\t$(ROOT)/$(MAPPER) -a $(MM2_OPT) $(MM2_IDX) $< 2> $@.log | gzip > $@});\n",
    "push(@mak, \"$pre.hap2.sam.gz:$hap[1]\");\n",
    "push(@mak, q{\t$(ROOT)/$(MAPPER) -a $(MM2_OPT) $(MM2_IDX) $< 2> $@.log | gzip > $@});\n",
    "push(@mak, \"\");\n",
    "\n",
    "push(@mak, \"$pre.hap1.bam:$pre.hap1.sam.gz\");\n",
    "push(@mak, q{\t$(ROOT)/k8 $(ROOT)/dipcall-aux.js samflt $< | $(ROOT)/samtools sort -m4G --threads 4 -o $@ -});\n",
    "push(@mak, \"$pre.hap2.bam:$pre.hap2.sam.gz\");\n",
    "push(@mak, q{\t$(ROOT)/k8 $(ROOT)/dipcall-aux.js samflt $< | $(ROOT)/samtools sort -m4G --threads 4 -o $@ -});\n",
    "push(@mak, \"\");\n",
    "\n",
    "push(@mak, \"$pre.pair.vcf.gz:$pre.hap1.bam $pre.hap2.bam\");\n",
    "push(@mak, q{\t$(ROOT)/htsbox pileup -q5 -evcf $(REF_FA) $^ | $(ROOT)/htsbox bgzip > $@});\n",
    "push(@mak, \"$pre.dip.vcf.gz:$pre.pair.vcf.gz\");\n",
    "\n",
    "my $pair_opt = defined($opts{a})? \"-a\" : \"\";\n",
    "if ($is_male) {\n",
    "\tpush(@mak, q{\t$(ROOT)/k8 $(ROOT)/dipcall-aux.js vcfpair} . qq{ $pair_opt -p $opts{x} } . q{$< | $(ROOT)/htsbox bgzip > $@});\n",
    "} else {\n",
    "\tpush(@mak, q{\t$(ROOT)/k8 $(ROOT)/dipcall-aux.js vcfpair} . qq{ $pair_opt } . q{$< | $(ROOT)/htsbox bgzip > $@});\n",
    "}\n",
    "push(@mak, \"\");\n",
    "\n",
    "push(@mak, \"$pre.hap1.var.gz:$pre.hap1.paf.gz\");\n",
    "push(@mak, q{\tgzip -dc $< | sort -k6,6 -k8,8n | $(ROOT)/k8 $(ROOT)/paftools.js call - 2> $@.vst | gzip > $@});\n",
    "push(@mak, \"$pre.hap2.var.gz:$pre.hap2.paf.gz\");\n",
    "push(@mak, q{\tgzip -dc $< | sort -k6,6 -k8,8n | $(ROOT)/k8 $(ROOT)/paftools.js call - 2> $@.vst | gzip > $@});\n",
    "push(@mak, \"$pre.hap1.bed:$pre.hap1.var.gz\");\n",
    "push(@mak, q{\tgzip -dc $< | grep ^R | cut -f2- > $@});\n",
    "push(@mak, \"$pre.hap2.bed:$pre.hap2.var.gz\");\n",
    "push(@mak, q{\tgzip -dc $< | grep ^R | cut -f2- > $@});\n",
    "push(@mak, \"$pre.dip.bed:$pre.hap1.bed $pre.hap2.bed\");\n",
    "if ($is_male) {\n",
    "\tmy $bedtk = '$(ROOT)/bedtk';\n",
    "\tmy @cmd;\n",
    "\tpush(@cmd, \"$bedtk isec $pre.hap1.bed $pre.hap2.bed | egrep -v '^(chr)?[XY]' > $pre.tmp.bed\"); # autosome\n",
    "\tpush(@cmd, \"$bedtk isec $pre.hap1.bed $pre.hap2.bed | $bedtk isec $opts{x} >> $pre.tmp.bed\");  # chrX PAR\n",
    "\tpush(@cmd, \"$bedtk sub $pre.hap2.bed $pre.hap1.bed | egrep '^(chr)?X' | $bedtk sub - $opts{x} >> $pre.tmp.bed\"); # chrX non-PAR\n",
    "\tpush(@cmd, \"$bedtk sub $pre.hap1.bed $pre.hap2.bed | egrep '^(chr)?Y' >> $pre.tmp.bed\");       # chrY non-PAR\n",
    "\tpush(@cmd, \"$bedtk sort $pre.tmp.bed > \" . '$@');\n",
    "\tpush(@mak, \"\\t\" . join(\"; \", @cmd));\n",
    "} else {\n",
    "\tpush(@mak, q{\t$(ROOT)/bedtk isec -m $^ > $@});\n",
    "}\n",
    "push(@mak, \"\");\n",
    "\n",
    "print(join(\"\\n\", @mak));\n",
    "\n",
    "exit(0);\n",
    "\n",
    "# auxiliary routines\n",
    "sub which {\n",
    "\tmy $file = shift;\n",
    "\tmy $path = (@_)? shift : $ENV{PATH};\n",
    "\treturn if (!defined($path));\n",
    "\tforeach my $x (split(\":\", $path)) {\n",
    "\t\t$x =~ s/\\/$//;\n",
    "\t\treturn \"$x/$file\" if (-x \"$x/$file\");\n",
    "\t}\n",
    "\treturn;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "f4afa7b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "!chmod u+x /results/output2/run_dipcall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1e86561",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "9091a476",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(\"/results/output2/\")\n",
    "for sample in sample_to_seqeunces:\n",
    "    os.system(f\"./run_dipcall -d grch38.mmidx DPC_{sample} ../grch38.fa {sample}_0_CMRG.fa {sample}_1_CMRG.fa > makefile_{sample}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "42058df2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Process is interrupted.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cd /results/output2/\n",
    "ls makefile_* | awk '{print \"make -f \"$1}' | parallel -j 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d3bf8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /results/\n",
    "tar zcvf all_CMRG_dipcalls_GRCH38.tgz output2\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}