4_merge_barcard_cell_barcodes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a57c7ba7-cd43-41cb-857a-23af0f5ce8a8",
   "metadata": {},
   "source": [
    "# 4. [OPTIONAL] Merging cell barcodes\n",
    "Based on the Otsu filtering in notebook #3, we now merge all the barcode sets with significant overlaps in the bams. This is necessary when performing analyses at .bam level (e.g. freemxulet) and you want to use the Barcard cell barcodes. Run the following code in the bash container found in `PUMATAC_dependencies`\n",
    "Open the bash container:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac3be90d-5c87-46a0-9add-73375d2cb2d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_path=PUMATAC_dependencies/cache/vibsinglecellnf-singlecelltoolkit-2022-07-07-0638c1d.img\n",
    "singularity exec --cleanenv -H $PWD -B /staging,/lustre1,/data,/vsc-hard-mounts,/scratch $img_path bash"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "597867e4-c338-4969-85ce-eff2c21453b9",
   "metadata": {},
   "source": [
    "In the container, run the following commands:  \n",
    "Generate fragments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f29d0d0-9c30-48a3-9fd8-82e1b3835971",
   "metadata": {},
   "outputs": [],
   "source": [
    "export POLARS_MAX_THREADS=8\n",
    "export OMP_NUM_THREADS=8\n",
    "\n",
    "for file in PUMATAC_tutorial_preprocessing_out/data/reports/barcard/*tomerge.tsv\n",
    "do\n",
    "    echo $file\n",
    "    sample=(`basename $file`)\n",
    "    sample=${sample%.tomerge.tsv}\n",
    "    # echo $sample\n",
    "    jaccardfile=$file\n",
    "    infile=PUMATAC_tutorial_preprocessing_out/data/fragments/$sample.fragments.raw.tsv.gz\n",
    "    outfile=PUMATAC_tutorial_preprocessing_out/data/fragments/$sample.fragments.tsv\n",
    "    outfile2=$outfile.gz\n",
    "\n",
    "\n",
    "    merge_cbs_over_jaccard_index_threshold.py -j $jaccardfile -i $infile -o $outfile &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b00ed5dc-ced7-4f7f-a841-952da8d034cd",
   "metadata": {},
   "source": [
    "Then bgzip the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "23d58377-7ee8-4f18-b195-0d7ca7c9662c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "full_preprocessing_out/data/fragments/OHS_s3atac_1.FULL.fragments.tsv\n",
      "exists\n",
      "full_preprocessing_out/data/fragments/OHS_s3atac_2.FULL.fragments.tsv\n",
      "exists\n"
     ]
    }
   ],
   "source": [
    "for file in PUMATAC_tutorial_preprocessing_out/data/fragments/*.fragments.tsv\n",
    "do\n",
    "    echo $file\n",
    "    bgzip -f -@ 8 $file &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17945f81-f5c1-486f-9e43-1f638bdf988e",
   "metadata": {},
   "source": [
    "And index the fragments files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da504ed8-2147-4e19-8b60-27fa4249c3b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for file in PUMATAC_tutorial_preprocessing_out/data/fragments/*.fragments.tsv.gz\n",
    "do\n",
    "    tabix -p bed $file &\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eb15bee-e3db-489a-98bf-4e87100ff202",
   "metadata": {},
   "source": [
    "If you want to process the `.bam` files further, add a new barcode tag `DB` to the corresponding .bam files, based on the barcode translate files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33a3ac70-5546-4aed-81fe-e9025edc11c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "prog=../0_resources/single_cell_toolkit_rust/target/release/add_new_bc_tag\n",
    "prog=/dodrio/scratch/projects/starting_2022_023/benchmark/scatac_benchmark/0_resources/single_cell_toolkit_rust/add_new_bc_tag_if_old_found\n",
    "for file in full_preprocessing_out/data/bam/UCS*possorted.bam\n",
    "do\n",
    "    echo $file\n",
    "    sample=(`basename $file`)\n",
    "    sample=${sample%_____R1.bwa.out.possorted.bam}\n",
    "    echo $sample\n",
    "    outfile=full_preprocessing_out/data/bam/$sample.possorted.dbfixed.bam\n",
    "    bc_map=full_preprocessing_out/data/fragments/$sample.fragments.tsv.original_CB_to_merged_CBs.tsv\n",
    "    ls $bc_map\n",
    "    if [ -f \"$outfile\" ]; then\n",
    "        echo \"$outfile exists, skipping\"\n",
    "    else \n",
    "        echo \"$outfile does not exist, generating!\"\n",
    "        $prog $file $outfile $bc_map CB DB &\n",
    "    fi\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c274cf3f-29b7-4d56-b58c-5fee9cac74d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "for file in full_preprocessing_out/data/bam/*.possorted.dbfixed.bam\n",
    "do\n",
    "    echo $file\n",
    "    outfile=$file.bai\n",
    "    if [ -f \"$outfile\" ]; then\n",
    "        echo \"$outfile exists, skipping\"\n",
    "    else \n",
    "        echo \"$outfile does not exist, generating!\"\n",
    "        samtools index -@ 8 $file &\n",
    "    fi\n",
    "done"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}