Unit2-module2-ChIPseq/ChIP-seq.practice4.6.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda\n",
    "mkdir ChIP-seq.practice4.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Use MACS2 with appropriate setting to call peaks for the file SRR3001750_srt_dupr.chr2.chr11.bam (downloaded from last exercise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO  @ Tue, 13 Nov 2018 09:50:08: \n",
      "# Command line: callpeak -t SRR3001750_srt_dupr.chr2.chr11.bam -f BAM -n SRR3001750_srt_dupr.chr2.chr11_MACS2 -g mm -p 1e-2 --broad --nomodel --bdg --shift 0 --extsize 147 --keep-dup all\n",
      "# ARGUMENTS LIST:\n",
      "# name = SRR3001750_srt_dupr.chr2.chr11_MACS2\n",
      "# format = BAM\n",
      "# ChIP-seq file = ['SRR3001750_srt_dupr.chr2.chr11.bam']\n",
      "# control file = None\n",
      "# effective genome size = 1.87e+09\n",
      "# band width = 300\n",
      "# model fold = [5, 50]\n",
      "# pvalue cutoff for narrow/strong regions = 1.00e-02\n",
      "# pvalue cutoff for broad/weak regions = 1.00e-01\n",
      "# qvalue will not be calculated and reported as -1 in the final output.\n",
      "# Larger dataset will be scaled towards smaller dataset.\n",
      "# Range for calculating regional lambda is: 10000 bps\n",
      "# Broad region calling is on\n",
      "# Paired-End mode is off\n",
      " \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:08: #1 read tag files... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:08: #1 read treatment tags... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:12:  1000000 \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #1 tag size is determined as 49 bps \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #1 tag size = 49 \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #1  total tags in treatment: 1675380 \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #1 finished! \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #2 Build Peak Model... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #2 Skipped... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #2 Use 147 as fragment length \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #3 Call peaks... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #3 Call broad peaks with given level1 -log10pvalue cutoff and level2: 2.00000, 1.00000... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:15: #3 Pre-compute pvalue-qvalue table... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:19: #3 In the peak calling step, the following will be performed simultaneously: \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:19: #3   Write bedGraph files for treatment pileup (after scaling if necessary)... SRR3001750_srt_dupr.chr2.chr11_MACS2_treat_pileup.bdg \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:19: #3   Write bedGraph files for control lambda (after scaling if necessary)... SRR3001750_srt_dupr.chr2.chr11_MACS2_control_lambda.bdg \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:19: #3 Call peaks for each chromosome... \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:41: #4 Write output xls file... SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.xls \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:42: #4 Write broad peak in broadPeak format file... SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.broadPeak \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:43: #4 Write broad peak in bed12/gappedPeak format file... SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.gappedPeak \n",
      "INFO  @ Tue, 13 Nov 2018 09:50:43: Done! \n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda/ChIP-seq.practice4.6\n",
    "macs2 callpeak -t SRR3001750_srt_dupr.chr2.chr11.bam -f BAM -n SRR3001750_srt_dupr.chr2.chr11_MACS2 -g mm -p 1e-2 --broad --nomodel --bdg --shift 0 --extsize 147 --keep-dup all "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Examine the out put and find out:\n",
    "* In broadPeaks file, how many peaks were identified in chr2? How many peaks were identified in chr11?\n",
    "* How many peaks in broadPeaks were identified with –log10pvalue  > 2?\n",
    "* Find the coordinates of peak with highest signal from treat_pileup.bdg file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr11\t3100344\t3104796\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_1\t35\t.\t2.80190\t6.25895\t3.59094\n",
      "chr11\t3106853\t3107002\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_2\t7\t.\t2.00168\t2.16457\t0.73649\n",
      "chr11\t3112273\t3113851\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_3\t6\t.\t2.14515\t1.99585\t0.69248\n",
      "chr11\t3114494\t3115651\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_4\t12\t.\t2.88502\t3.18291\t1.21582\n",
      "chr11\t3117596\t3117743\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_5\t5\t.\t1.94322\t1.66785\t0.55325\n",
      "chr11\t3118845\t3119359\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_6\t6\t.\t1.98141\t1.83885\t0.62141\n",
      "chr11\t3125426\t3129535\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_7\t22\t.\t1.75859\t4.70931\t2.23259\n",
      "chr11\t3130943\t3135775\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_8\t9\t.\t1.41238\t2.53377\t0.96321\n",
      "chr11\t3136482\t3137703\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_9\t8\t.\t1.43226\t2.48031\t0.88682\n",
      "chr11\t3139576\t3140021\tSRR3001750_srt_dupr.chr2.chr11_MACS2_peak_10\t4\t.\t1.28025\t1.45590\t0.48721\n",
      "   73346\n",
      "   47241\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda/ChIP-seq.practice4.6\n",
    "head SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.broadPeak\n",
    "awk -F \"\\t\" '{if ($1==\"chr2\") {print $0}}' SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.broadPeak | wc -l\n",
    "awk -F \"\\t\" '{if ($1==\"chr11\") {print $0}}' SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.broadPeak | wc -l\n",
    "# Can you use grep here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   21329\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda/ChIP-seq.practice4.6\n",
    "awk -F \"\\t\" '{if ($8>2) {print $0}}' SRR3001750_srt_dupr.chr2.chr11_MACS2_peaks.broadPeak | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr11\t0\t3099930\t0.00000\n",
      "chr11\t3099930\t3099938\t1.00000\n",
      "chr11\t3099938\t3100077\t2.00000\n",
      "chr11\t3100077\t3100085\t1.00000\n",
      "chr11\t3100085\t3100152\t0.00000\n",
      "chr11\t3100152\t3100177\t1.00000\n",
      "chr11\t3100177\t3100184\t2.00000\n",
      "chr11\t3100184\t3100254\t3.00000\n",
      "chr11\t3100254\t3100287\t4.00000\n",
      "chr11\t3100287\t3100299\t5.00000\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda/ChIP-seq.practice4.6\n",
    "head SRR3001750_srt_dupr.chr2.chr11_MACS2_treat_pileup.bdg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda/ChIP-seq.practice4.6\n",
    "sort -nrk4,4 SRR3001750_srt_dupr.chr2.chr11_MACS2_treat_pileup.bdg > SRR3001750_srt_dupr.chr2.chr11_MACS2_treat_pileup_ranked.bdg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr2\t98667174\t98667189\t294.00000\n",
      "chr2\t98666849\t98666892\t294.00000\n",
      "chr2\t98666686\t98666702\t294.00000\n",
      "chr2\t98667189\t98667190\t293.00000\n",
      "chr2\t98667172\t98667174\t293.00000\n",
      "chr2\t98666892\t98666895\t293.00000\n",
      "chr2\t98666702\t98666849\t293.00000\n",
      "chr2\t98666683\t98666686\t293.00000\n",
      "chr2\t98667190\t98667191\t292.00000\n",
      "chr2\t98667171\t98667172\t292.00000\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cd /Users/yolandatiao/Documents/0_Bioinformatics2017/2018_Bioinformatics/Applied-Bioinformatics-HW-Yolanda/ChIP-seq.practice4.6\n",
    "head SRR3001750_srt_dupr.chr2.chr11_MACS2_treat_pileup_ranked.bdg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}