{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# __DAY3 - MAPPING PRACTICE__" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this practice, we are going to map short reads against a reference. To know, how well do the reads align back to the reference, we use bwa-mem2 and samtools to assess the basic alignment statistics.\n", "\n", "In this exercise, we will use :\n", "* the ___Bathycoccus prasinos_ genome__ - File : `GCA_002220235.1_ASM222023v1_genomic.fna`\n", "* the __ILLUMINA READS of the 1613 sample__ - Files `1613_R1.fastq.gz`, `1613_R2.fastq.gz`\n", "\n", "The tool bwa needs 2 steps: \n", "- **Reference indexing**: `bwa-mem2 index reference`\n", "- **Mapping in itself**: `bwa-mem2 mem -R READGROUP [options] reference fastq1 fastq2 > out.sam`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reference indexation - `bwa_mem2 index`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before mapping, we need index reference file! \n", "\n", "#### Go into REF directory and list the content of the directory" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/home/jovyan/SV_DATA/REF\n", "total 561996\n", "-rw-r--r-- 1 jovyan users 380710511 Oct 31 2011 all.chrs.con\n", "-rw-r--r-- 1 jovyan users 81498659 Feb 7 2012 all.gff3\n", "-rw------- 1 jovyan users 15263945 Sep 23 01:29 GCA_002220235.1_ASM222023v1_genomic.fna\n", "-rw------- 1 jovyan users 5886829 Sep 23 15:25 GCA_002220235.1_ASM222023v1_genomic.gff\n", "-rw-r--r-- 1 jovyan users 3768582 Sep 24 10:45 GCA_002220235.1_ASM222023v1_genomic.fna.pac\n", "-rw-r--r-- 1 jovyan users 1587 Sep 24 10:45 GCA_002220235.1_ASM222023v1_genomic.fna.ann\n", "-rw-r--r-- 1 jovyan users 308 Sep 24 10:45 GCA_002220235.1_ASM222023v1_genomic.fna.amb\n", "-rw-r--r-- 1 jovyan users 30148640 Sep 24 10:45 GCA_002220235.1_ASM222023v1_genomic.fna.0123\n", "-rw-r--r-- 1 jovyan users 48991633 Sep 24 10:45 GCA_002220235.1_ASM222023v1_genomic.fna.bwt.2bit.64\n", "-rw-r--r-- 1 jovyan users 681 Sep 25 19:37 GCA_002220235.1_ASM222023v1_genomic.fna.fai\n", "-rw-r--r-- 1 jovyan users 65 Oct 2 14:17 qtl.bed\n", "-rw-r--r-- 1 jovyan users 594953 Oct 3 13:16 qtl.fasta\n", "-rw-r--r-- 1 jovyan users 604866 Oct 3 13:18 qtl.formatseqtk.fasta\n", "-rw-r--r-- 1 jovyan users 353 Oct 3 13:19 all.chrs.con.fai\n", "-rw-r--r-- 1 jovyan users 7957814 Oct 3 13:23 all_genes.gff3\n", "-rw-r--r-- 1 jovyan users 14204 Oct 3 13:27 qtl.intesect.genes.bed\n" ] } ], "source": [ "cd /home/jovyan/SV_DATA/REF/\n", "pwd\n", "ls -lrt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Print the first lines of the file that contains the genome sequence (fasta format)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">FO082278.1 Bathycoccus prasinos genomic : Chromosome_1\n", "TTGCATAGGTTTCTAGGGTCTGGGATTTTTGGGATTTTCCGAAAACTCCTTTCTCTAAACTGGACGgtctcggcggcgac\n", "gagacTCCCGGGTTGGGAATGTTCTTTTGCGCTGCAATCATATCCTGCTCCGGATACCGAATCGTGCGCTGCGAAATCGT\n", "CTCCGCACTTGAAGGCTGTTTCTCACCATTACTATTATTAACGCTATCTCGAACTTCCTGCTTCGGTGCCGCCGGACCTA\n", "CTCTGTTATTTCCCACGAGCTGTTgtttctcgctctctttgcTCACTTTTACTACATCACTCGTCGCGTCGAGAACATGC\n" ] } ], "source": [ "head /home/jovyan/SV_DATA/REF/GCA_002220235.1_ASM222023v1_genomic.fna -n 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Before mapping we need index reference file! Check bwa-mem2 index command line." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Usage: bwa-mem2 \n", "Commands:\n", " index create index\n", " mem alignment\n", " version print version number\n" ] }, { "ename": "", "evalue": "1", "output_type": "error", "traceback": [] } ], "source": [ "bwa-mem2" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Usage: bwa index [-p prefix] \n" ] }, { "ename": "", "evalue": "1", "output_type": "error", "traceback": [] } ], "source": [ "#PUT COMMAND TO GET HEP FOR bwa-mem2 index\n", "bwa-mem2 index" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND TO CREATE INDEX" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bwa_index] Pack FASTA... 0.06 sec\n", "init ticks = 1077859628\n", "ref seq len = 30148640\n", "binary seq ticks = 553871378\n", "build index ticks = 12798787546\n", "ref_seq_len = 30148640\n", "count = 0, 7834025, 15074320, 22314615, 30148640\n", "BWT[28778507] = 4\n", "CP_SHIFT = 5, CP_MASK = 31\n", "sizeof CP_OCC = 64\n", "max_occ_ind = 942145\n", "ref_seq_len = 30148640\n", "count = 0, 7834025, 15074320, 22314615, 30148640\n", "BWT[28778507] = 4\n", "CP_SHIFT = 6, CP_MASK = 63\n", "sizeof CP_OCC = 64\n", "max_occ_ind = 471072\n" ] } ], "source": [ "bwa-mem2 index GCA_002220235.1_ASM222023v1_genomic.fna" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Check that the indexes have been created" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 896904\n", "-rw-r--r-- 1 jovyan users 380710511 Oct 31 2011 all.chrs.con\n", "-rw-r--r-- 1 jovyan users 81498659 Feb 7 2012 all.gff3\n", "-rw------- 1 jovyan users 15263945 Sep 23 01:29 GCA_002220235.1_ASM222023v1_genomic.fna\n", "-rw------- 1 jovyan users 5886829 Sep 23 15:25 GCA_002220235.1_ASM222023v1_genomic.gff\n", "-rw-r--r-- 1 jovyan users 681 Sep 25 19:37 GCA_002220235.1_ASM222023v1_genomic.fna.fai\n", "-rw-r--r-- 1 jovyan users 65 Oct 2 14:17 qtl.bed\n", "-rw-r--r-- 1 jovyan users 594953 Oct 3 13:16 qtl.fasta\n", "-rw-r--r-- 1 jovyan users 604866 Oct 3 13:18 qtl.formatseqtk.fasta\n", "-rw-r--r-- 1 jovyan users 353 Oct 3 13:19 all.chrs.con.fai\n", "-rw-r--r-- 1 jovyan users 7957814 Oct 3 13:23 all_genes.gff3\n", "-rw-r--r-- 1 jovyan users 14204 Oct 3 13:27 qtl.intesect.genes.bed\n", "-rw-r--r-- 1 jovyan users 3768582 Oct 3 13:47 GCA_002220235.1_ASM222023v1_genomic.fna.pac\n", "-rw-r--r-- 1 jovyan users 1587 Oct 3 13:47 GCA_002220235.1_ASM222023v1_genomic.fna.ann\n", "-rw-r--r-- 1 jovyan users 308 Oct 3 13:47 GCA_002220235.1_ASM222023v1_genomic.fna.amb\n", "-rw-r--r-- 1 jovyan users 30148640 Oct 3 13:47 GCA_002220235.1_ASM222023v1_genomic.fna.0123\n", "-rw-r--r-- 1 jovyan users 211040597 Oct 3 13:47 GCA_002220235.1_ASM222023v1_genomic.fna.bwt.8bit.32\n", "-rw-r--r-- 1 jovyan users 180891925 Oct 3 13:47 GCA_002220235.1_ASM222023v1_genomic.fna.bwt.2bit.64\n" ] } ], "source": [ "ls -lrt" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Run the mapping with `bwa mem` for only one sample (sample 1613)\n", "\n", "\n", "#### Go into the directory SV_DATA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/home/jovyan/SV_DATA/results\n", "/home/jovyan/SV_DATA\n" ] } ], "source": [ "pwd\n", "cd /home/jovyan/SV_DATA\n", "pwd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create a direcctory `RESULTS` and go into the directory " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "### UPDATING\n", "mkdir -p RESULTS" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/home/jovyan/SV_DATA/RESULTS\n" ] } ], "source": [ "cd RESULTS\n", "pwd" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#PUT COMMAND TO GET HEP FOR bwa-mem2 mem" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-----------------------------\n", "Executing in AVX2 mode!!\n", "-----------------------------\n", "Usage: bwa2 mem [options] [in2.fq]\n", "Options:\n", " Algorithm options:\n", " -o STR Output SAM file name\n", " -t INT number of threads [1]\n", " -k INT minimum seed length [19]\n", " -w INT band width for banded alignment [100]\n", " -d INT off-diagonal X-dropoff [100]\n", " -r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [1.5]\n", " -y INT seed occurrence for the 3rd round seeding [20]\n", " -c INT skip seeds with more than INT occurrences [500]\n", " -D FLOAT drop chains shorter than FLOAT fraction of the longest overlapping chain [0.50]\n", " -W INT discard a chain if seeded bases shorter than INT [0]\n", " -m INT perform at most INT rounds of mate rescues for each read [50]\n", " -S skip mate rescue\n", " -o output file name missing\n", " -P skip pairing; mate rescue performed unless -S also in use\n", "Scoring options:\n", " -A INT score for a sequence match, which scales options -TdBOELU unless overridden [1]\n", " -B INT penalty for a mismatch [4]\n", " -O INT[,INT] gap open penalties for deletions and insertions [6,6]\n", " -E INT[,INT] gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1]\n", " -L INT[,INT] penalty for 5'- and 3'-end clipping [5,5]\n", " -U INT penalty for an unpaired read pair [17]\n", "Input/output options:\n", " -p smart pairing (ignoring in2.fq)\n", " -R STR read group header line such as '@RG\\tID:foo\\tSM:bar' [null]\n", " -H STR/FILE insert STR to header if it starts with @; or insert lines in FILE [null]\n", " -j treat ALT contigs as part of the primary assembly (i.e. ignore .alt file)\n", " -v INT verbose level: 1=error, 2=warning, 3=message, 4+=debugging [3]\n", " -T INT minimum score to output [30]\n", " -h INT[,INT] if there are 80% of the max score, output all in XA [5,200]\n", " -a output all alignments for SE or unpaired PE\n", " -C append FASTA/FASTQ comment to SAM output\n", " -V output the reference FASTA header in the XR tag\n", " -Y use soft clipping for supplementary alignments\n", " -M mark shorter split hits as secondary\n", " -I FLOAT[,FLOAT[,INT[,INT]]]\n", " specify the mean, standard deviation (10% of the mean if absent), max\n", " (4 sigma from the mean if absent) and min of the insert size distribution.\n", " FR orientation only. [inferred]\n", "Note: Please read the man page for detailed description of the command line and options.\n", "No. of OMP threads: 0\n", "Processor is runnig @2295.132606 MHz\n", "Runtime profile:\n", "\n", "\t Time taken for main_mem function: 0.00 Sec\n", "\n", "\tIO times (sec) :\n", "\tReading IO time (reads) avg: 0.00, (0.00, 0.00)\n", "\tWriting IO time (SAM) avg: 0.00, (0.00, 0.00)\n", "\tReading IO time (Reference Genome) avg: 0.00, (0.00, 0.00)\n", "\tIndex read time avg: 0.00, (0.00, 0.00)\n", "\n", "\tOverall time (sec) (Excluding Index reading time):\n", "\tPROCESS() (Total compute time + (read + SAM) IO time) : 0.00\n", "\tMEM_PROCESS_SEQ() (Total compute time (Kernel + SAM)), avg: 0.00, (0.00, 0.00)\n", "\n", "\t SAM Processing time (sec):\n", "\t--WORKER_SAM avg: 0.00, (0.00, 0.00)\n", "\n", "\tKernels' compute time (sec):\n", "\tTotal kernel (smem+sal+bsw) time avg: 0.00, (0.00, 0.00)\n", "\t\tSMEM compute avg: -nan, (0.00, 435704.67)\n", "\t\tSAL compute avg: -nan, (0.00, 435704.67)\n", "\t\tBSW time, avg: -nan, (0.00, 435704.67)\n", "\n", "\tTotal allocs: 0 = 0 out total requests: 0, Rate: -nan\n", "\n", "Important parameter settings: \n", "\tBATCH_SIZE: 512\n", "\tMAX_SEQ_LEN_REF: 256\n", "\tMAX_SEQ_LEN_QER: 128\n", "\tMAX_SEQ_LEN8: 128\n", "\tSEEDS_PER_READ: 500\n", "\tSIMD_WIDTH8 X: 32\n", "\tSIMD_WIDTH16 X: 16\n", "\tAVG_SEEDS_PER_READ: 64\n" ] } ], "source": [ "bwa-mem2 mem" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-----------------------------\n", "Executing in AVX2 mode!!\n", "-----------------------------\n", "Ref file: ../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "Entering FMI_search\n", "reference seq len = 30148641\n", "count\n", "0,\t1\n", "1,\t7834026\n", "2,\t15074321\n", "3,\t22314616\n", "4,\t30148641\n", "\n", "Reading other elements of the index from files ../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "prefix: ../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "[M::bwa_idx_load_ele] read 0 ALT contigs\n", "Done reading Index!!\n", "Reading reference genome..\n", "Binary seq file = ../REF/GCA_002220235.1_ASM222023v1_genomic.fna.0123\n", "Reference genome size: 30148640 bp\n", "Done readng reference genome !!\n", "\n", "[0000] 1: Calling process()\n", "\n", "Threads used (compute): 8\n", "Info: projected #read in a task: 529811\n", "------------------------------------------\n", "Memory pre-allocation for chaining: 1114.7223 MB\n", "Memory pre-allocation for BSW: 1916.9362 MB\n", "Memory pre-allocation for BWT: 618.5134 MB\n", "------------------------------------------\n", "No. of pipeline threads: 2\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000217, nseq: 320478\n", "\t[0000][ M::kt_pipeline] read 320478 sequences (80000217 bp)...\n", "[0000] 2. Calling mem_process_seqs.., task: 0\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000178, nseq: 320478\n", "\t[0000][ M::kt_pipeline] read 320478 sequences (80000178 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320478\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (6, 146360, 9, 7)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (228, 278, 344)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 576)\n", "[M::mem_pestat] mean and std.dev: (290.09, 87.34)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 692)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320478 reads in 47.920 CPU sec, 6.018 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 1\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000142, nseq: 320476\n", "\t[0000][ M::kt_pipeline] read 320476 sequences (80000142 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320478\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (7, 146471, 14, 7)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (227, 277, 343)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 575)\n", "[M::mem_pestat] mean and std.dev: (289.75, 87.63)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 691)\n", "[M::mem_pestat] analyzing insert size distribution for orientation RF...\n", "[M::mem_pestat] (25, 50, 75) percentile: (60, 334, 2198)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 6474)\n", "[M::mem_pestat] mean and std.dev: (1083.00, 1224.54)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 8612)\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[M::mem_pestat] skip orientation RF\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320478 reads in 45.630 CPU sec, 5.652 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 2\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000457, nseq: 320478\n", "\t[0000][ M::kt_pipeline] read 320478 sequences (80000457 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320476\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (13, 146615, 11, 9)\n", "[M::mem_pestat] analyzing insert size distribution for orientation FF...\n", "[M::mem_pestat] (25, 50, 75) percentile: (3706, 4374, 7453)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 14947)\n", "[M::mem_pestat] mean and std.dev: (4598.31, 2463.35)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 18694)\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (227, 277, 343)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 575)\n", "[M::mem_pestat] mean and std.dev: (289.61, 87.35)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 691)\n", "[M::mem_pestat] analyzing insert size distribution for orientation RF...\n", "[M::mem_pestat] (25, 50, 75) percentile: (156, 217, 2198)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 6282)\n", "[M::mem_pestat] mean and std.dev: (1123.36, 1223.29)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 8324)\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[M::mem_pestat] skip orientation FF\n", "[M::mem_pestat] skip orientation RF\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320476 reads in 46.026 CPU sec, 5.680 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 3\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 47353187, nseq: 189696\n", "\t[0000][ M::kt_pipeline] read 189696 sequences (47353187 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320478\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (6, 146206, 5, 5)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (227, 276, 342)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 572)\n", "[M::mem_pestat] mean and std.dev: (288.48, 86.88)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 687)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320478 reads in 45.713 CPU sec, 5.704 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 4\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 0, nseq: 0\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 189696\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 86576, 7, 5)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (227, 277, 343)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 575)\n", "[M::mem_pestat] mean and std.dev: (289.32, 87.35)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 691)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 189696 reads in 27.671 CPU sec, 3.486 real sec\n", "[0000] read_chunk: 80000000, work_chunk_size: 0, nseq: 0\n", "[0000] Computation ends..\n", "No. of OMP threads: 8\n", "Processor is runnig @2295.048362 MHz\n", "Runtime profile:\n", "\n", "\t Time taken for main_mem function: 28.08 Sec\n", "\n", "\tIO times (sec) :\n", "\tReading IO time (reads) avg: 4.88, (4.88, 4.88)\n", "\tWriting IO time (SAM) avg: 1.94, (1.94, 1.94)\n", "\tReading IO time (Reference Genome) avg: 0.02, (0.02, 0.02)\n", "\tIndex read time avg: 0.18, (0.18, 0.18)\n", "\n", "\tOverall time (sec) (Excluding Index reading time):\n", "\tPROCESS() (Total compute time + (read + SAM) IO time) : 27.88\n", "\tMEM_PROCESS_SEQ() (Total compute time (Kernel + SAM)), avg: 26.53, (26.53, 26.53)\n", "\n", "\t SAM Processing time (sec):\n", "\t--WORKER_SAM avg: 3.84, (3.84, 3.84)\n", "\n", "\tKernels' compute time (sec):\n", "\tTotal kernel (smem+sal+bsw) time avg: 22.38, (22.38, 22.38)\n", "\t\tSMEM compute avg: 8.15, (8.22, 8.12)\n", "\t\tSAL compute avg: 0.85, (0.88, 0.81)\n", "\t\tBSW time, avg: 12.87, (12.93, 12.77)\n", "\n", "\tTotal allocs: 1941303 = 1941303 out total requests: 2368737, Rate: 0.82\n", "\n", "Important parameter settings: \n", "\tBATCH_SIZE: 512\n", "\tMAX_SEQ_LEN_REF: 256\n", "\tMAX_SEQ_LEN_QER: 128\n", "\tMAX_SEQ_LEN8: 128\n", "\tSEEDS_PER_READ: 500\n", "\tSIMD_WIDTH8 X: 32\n", "\tSIMD_WIDTH16 X: 16\n", "\tAVG_SEEDS_PER_READ: 64\n" ] } ], "source": [ "bwa-mem2 mem -R \"@RG\\tID:1613\\tSM:1613\" -M -t 8 ../REF/GCA_002220235.1_ASM222023v1_genomic.fna ../SHORT_READS/1613_R1.fastq.gz ../SHORT_READS/1613_R2.fastq.gz > 1613.sam" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "#### Check that the file `.sam` have been created by `bwa mem` " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 902M\n", "-rw-r--r-- 1 jovyan users 902M Oct 3 14:01 1613.sam\n" ] } ], "source": [ "ls -lh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Display the first and the end of the sam file just created" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMANDS" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "@SQ\tSN:FO082278.1\tLN:1352724\n", "@SQ\tSN:FO082277.1\tLN:1122692\n", "@SQ\tSN:FO082276.1\tLN:1091008\n", "@SQ\tSN:FO082275.1\tLN:1037991\n", "@SQ\tSN:FO082274.1\tLN:1019276\n", "@SQ\tSN:FO082273.1\tLN:989707\n", "@SQ\tSN:FO082272.1\tLN:955652\n", "@SQ\tSN:FO082271.1\tLN:937610\n", "@SQ\tSN:FO082270.1\tLN:895536\n", "@SQ\tSN:FO082269.1\tLN:794368\n", "@SQ\tSN:FO082268.1\tLN:741603\n", "@SQ\tSN:FO082267.1\tLN:712459\n", "@SQ\tSN:FO082266.1\tLN:708035\n", "@SQ\tSN:FO082265.1\tLN:663424\n", "@SQ\tSN:FO082264.1\tLN:519835\n", "@SQ\tSN:FO082263.1\tLN:494108\n", "@SQ\tSN:FO082262.1\tLN:465570\n", "@SQ\tSN:FO082261.1\tLN:310170\n", "@SQ\tSN:FO082260.1\tLN:146238\n", "@SQ\tSN:FO082258.2\tLN:43614\n", "@SQ\tSN:FO082259.2\tLN:72700\n", "@RG\tID:1613\tSM:1613\n", "@PG\tID:bwa\tPN:bwa\tVN:2.0pre2\tCL:bwa-mem2 mem -R @RG\\tID:1613\\tSM:1613 -M -t 8 ../REF/GCA_002220235.1_ASM222023v1_genomic.fna ../SHORT_READS/1613_R1.fastq.gz ../SHORT_READS/1613_R2.fastq.gz\n", "M02713:279:GW200319000:1:2106:15454:1334\t83\tFO082272.1\t908208\t60\t18S232M\t=\t908208\t-232\tCGTGTGCTCTTCCGATCTAACCGAGCGTGCCTTTCTTCTTTTTTCTGTAAAAGGAAGAACTAGTAGTAGTCGTAGCGTCGCCGTCTTCGTCTTCCATCGTCAAAGCTTCTTTCCTGCGTCTCTCCACGGTGCTCGCTCGAACCCCACCGACGGTTATCCCTTGATTATCGCTCGTTTTCCGTCTCCCTCGTTTTTGTTGCACGACGGTCCATCCGTCGTCTTCCTCTTCTCTCGCGCGTTTCCGTTTTTC\t-FFFF..EFFFDEFFFFDFFB;FFFFFFEA/FFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFDA?FFFFFFFFFFFFFFFFFBFDFFFFFFFFFFFFFBGEGGAGGGGGEGGEDCE;AEHGD?F@C/G?HCGCFDEFHDGHG>GHGH1CGHFHHEFFF>FE>/EHGHEHHGGEHGDDA@FH/@FGGGGGCC?@@C@FECFHBGBH/DHFGHHE?AGCAE.CH0.9../9/BFGC.=EA9CFFBF0;.@9@AB@BFBAB-;@-@9A9/9FF9/B/BFFF@@@=9@BBABFFEAFBF//9-9A---;B-AB>@\tNM:i:2\tMD:Z:71C125T34\tAS:i:222\tXS:i:19\tRG:Z:1613\n", "M02713:279:GW200319000:1:2106:15870:1335\t99\tFO082278.1\t612747\t60\t250M\t=\t612783\t286\tGTCTGGATCGAGTGCCTTCGGTAAACGGCGGGAATTTCGAGGGCGTGGACATCCCATCGATTCAAACAGTCGCCGAGCGATGGGCATCGAGAAACGGCTGCCAAACGGAAGAATCGGTGGTTGTTGATGGAGCGCTGAATTTGCGCACGAAGCTTTTCGCGAATGCGCCGGATACTTCCTCAGAAACGTGGACGAATGGTACGTGCCAAGAAGGAACAGTTGTTGAGAATTGGGTAATCGACGGGATGAA\tA>AABFFDC?AA44AGFGGGGFHDF5AAEEGGGAAGHFFGHGCDEE/EAG3B@?G3@31?/GEFGHHHB4FEEECGCDCD/B?GHHC0/AA.DFFFF/?FBFA9FBB.DFE.DFB;.FFGDFGFFFGBADFGE?EEFFFGF;/.DFGA0CC0HFD?CCDCC<@DGG0HGFD@\tNM:i:4\tMD:Z:13T39G109C62A23\tAS:i:230\tXS:i:24\tRG:Z:1613\n", "M02713:279:GW200319000:1:2106:14993:1368\t83\tFO082272.1\t365940\t43\t98S144M1D8M\t=\t365940\t-153\tAATGTGATAGAGTGAGCTTGGTTTTTTTTTTTCAAGCAGAAGACGGCATACGAGATCGTAGGTTGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTCATATGTTTAAGTCCATATGCGTATTTTAGTCGGCGGCACTACTAAAACCCGTTTCAAAGTCTCAGTCCTGCTGAGGCGAAGTCTTATGTTTCCATGAGACGGCGTGGTGGACCGTCGCGACTCCGTTCCGCGGGGCTCGAAACCGCTTTTT\t///;/;//////9//......-FDCFFFFFBFFFFFFFFFFFFFFFFFFFFFA.FFFFFFFFFFFFFFFFFFFFFFFFFFFGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGHGGGGHHHHHHHGGGGGGGGHHHHGHHHHGEF>GGGHHHHHHHHHDGHHHHGHHHHHHGGGGGHHHHHHHHHHHFGHGHHHHHGGGGGGHHHHHGGGGGGGGGGHGGGGGGGGGGGGGGGGGGGGBBBBBFBBBBBB\tNM:i:1\tMD:Z:144^A8\tAS:i:145\tXS:i:140\tRG:Z:1613\n" ] } ], "source": [ "head 1613.sam -n 30" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "M02713:279:GW200319000:1:1110:15626:29256\t99\tFO082262.1\t55546\t60\t156M94S\t=\t55546\t156\tCGTCGATAGATCCATGTTGAAATTCATTGCATAAGAAAACATGCGCGTCATATTTTTCACGTTGGATACGTCCCAACGATCCAAAAAGCCATTGAAGCTGGAACGCGATTCAAACGTATTATTCATCTGCGTCACTTGGCTGACGTCCCAGTACGGAGATCGGAAGAGCACACGTCTGAACTCCAGTCACAACCTGCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAGCGTGACTAAACACGG\tAAAABDABBFFFGGGGGGGGGGHHHHHHHHHHHHHHFGHGHHHHEAEEFGHHHHHHHFGHFGHHHHHHFFHGCEGHGGFHHGHGHHHDFHHHHHHHHHHGHGHHEEFGGHHHHHHGGHHHHHFGGGHHHHGGGGGHFHHHHGHHGGGGGHHGHHGGDD@HHFFGGGHGCGHHHGHHHGHHHHHFGGGHFHHHGGHGGGFA.FFGADFFGBBGGGGGGGFGFBFF/BFFAF-@9-......9//99;/;A-\tNM:i:2\tMD:Z:34T71C49\tAS:i:146\tXS:i:133\tRG:Z:1613\tXA:Z:FO082262.1,+58333,153M97S,4;\n", "M02713:279:GW200319000:1:1110:15626:29256\t147\tFO082262.1\t55546\t60\t94S156M\t=\t55546\t-156\tTGTATCTCAGGTTTTTTTTTTCTAACTCAGACGGACACCGCCCATATCTACACAAGCGGAAACACTCTTTCCCTACACGACGCTCTTCCGATCTCGTCGATAGATCCATGTTGAAATTCATTGCATAAGAAAACATGCGCGTCATATTTTTCACGTTGGATACGTCCCAACGATCCAAAAAGCCATTGAAGCTGGAACGCGATTCAAACGTATTATTCATCTGCGTCACTTGGCTGACGTCCCAGTACGG\t/;///;9./..-9@--:.9/9/9////.....//9--9-..;09090;090;A99CC0B/C0FGFFFEGF;0C:.DGGE;GFHGGGCGFHEAGCGCEGHHHHHDGFHGHHHHHHHHHHHHHFGHHFHHHHHHHGGGGGEGFHHHGHHHHGGHHGHHHEHGHGGHHGHHHHGHGHGGHHHHHHHHHHHHHHHHHGFAEEGHGHGHGHFHHHHHHHHHFHHGGGGGHGHHGGGFGEEGGF?FBFFA3AAA?A\tNM:i:2\tMD:Z:34T71C49\tAS:i:146\tXS:i:133\tRG:Z:1613\tXA:Z:FO082262.1,-58333,94S156M,5;\n", "M02713:279:GW200319000:1:1110:16510:29256\t83\tFO082273.1\t820998\t60\t36S213M\t=\t820998\t-213\tTTGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTAGTCCACGAAAGGGATGTTTGAAGTTTCCGATCACGATTTGCTGTTTTCGCAAGGTGACTCGACGGCAGCGGGTAAAGGAGCAGCTAAAGTACCGGATTCAAAGTTAGGCGGCGGGAAGAAGAAAGGCTCGTCGAAGACCGGTAGCGGAGGAAGTGGCGGTGGCGGCATTCTCGGTAGTAAAGGAAACGCGAACAACGCGAATTGGGGCGGCG\tBBFFBFFF;B/FFFFFFFFFFFFFFF??DFEFF/FFFFFDEAFFFFFFFFFFFFFFFFFFFFFFFFB/FFFFFFFFFFFFF>@-GGGGGGFGGGGBCGGAGGFGFAFHHGHHHHGHHFHHFHHHHG@A?GHHHHHHHHHHHGGGGGGGGHHHHHHHHHHHGGGGGGGHHHHGGGCGHGGGGGHHHHHHHGGGGFHGGFGGHHGGE?CGFHGGHHHGHHGHFE?GGGGFFCGEGGGGGFBBBB?BBAAAA\tNM:i:1\tMD:Z:65A147\tAS:i:208\tXS:i:0\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:16510:29256\t163\tFO082273.1\t820998\t60\t191M58S\t=\t820998\t213\tAGTCCACGAAAGGGATGTTTGAAGTTTCCGATCACGATTTGCTGTTTTCGCAAGGTGACTCGACGGCAGCGGGTAAAGGAGCAGCTAAAGTACCGGATTCAAAGTTAGGCGGCGGGAAGAAGAAAGGCTCGTCGAAGACCGGTAGCGGAGGAAGTGGCGGTGGGGGCGGTCTCTGTAGTAAAAGAAACGCGGAAAAAACGAATTGGGGGGGGGAGAGTGGCATAGAGTAGGGGGGGGAGAGAAGGATTT\t3AB?AFFABAABCF4GGGFGGGHGGHFGGDGGCFFDDHHGHHFHHHHHGGEGGGGEGHGFGGHGGGGGGGGGGEEHFHGHHGFFHHFBGEDGHFGGGFFDHHGHDGHHEFDCCGGGGGGHHHHCGBGEFHGGGG@EFGGGGB-ADE@;B-@DFFFFB.B=-9.--9---;-.////////;:/9..9/.9-----..-.;..;A/9A-----9--..//...////99/99..9--9--;..//..///\tNM:i:6\tMD:Z:65A97C3A0T4G8G8\tAS:i:161\tXS:i:21\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:15799:29276\t99\tFO082271.1\t640207\t60\t10M6D110M3D129M\t=\t640396\t439\tCAGTTTCTTGAACAACAACAACAACGATGCTTTTCTCGCTGCTCTTATTAAAGTCGTCGTGATACTAATAATGGGAGGAGATGAAAAAGACGATTTGAAAAAGTACCACCACCGTCGTCATCGTCGAGGAGGTGAGACGACTTTTAGAAAGCGATTCGCGCTCTACTTCGTGTGGTTCCTCGCGCTTTCCGTGTGGACGAGTTGGCAGATTTCCAATCGGTTCATGACGAAGATGACGATGCACGAGGA\t>AA>AFFFFFFFGGGAEGGGG1FH0EFHG?GGFHHHFF0FEEHHHHHHHHHHHHFCFFEEGHGHHHHHGHHHHHHHGGGGHHHHHHHHGGFEGEHGHHHHFHGHHFFHEEFCFEGG@EFGHHHHFCGGGGGGGHGFFFGG?FGHGHHGHHHGGGHGHGGGC?FEHHHHHGCECGGHE;GGFGGG@GGGGGGFEFGBFG@@@FFFEFFBFFFFFFFFFFFFFFFFFFFFFE@BFFFFFFFFFFFEEF@@@\tNM:i:11\tMD:Z:10^AACAAC49C20G39^TCC129\tAS:i:220\tXS:i:19\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:15799:29276\t147\tFO082271.1\t640396\t60\t250M\t=\t640207\t-439\tCGCGCTTTCCGTGTGGACGAGTTGGCAGATTTCCAATCGGTTCATGACGAAGATGACGATGCACGAGGACGAGAGCAACACCTTTCAAAACAACAACAAGATCGTCCTCAGTCACGAACGATGCGACGACTTCGATTTGGACGAGAACTTTTCTTTTTCGACGTCTTCCCCGACGAAATGCACGTGCGAAAGAGGCTACGGAGGTGCCAAATGCCGAGACGCGATGTGCGCGCGCCAGTGCGCCCCGCAC\tB;A;9-?B;?DFFD.DFD9;.:;/FFFFFFFFF.A=@;.9FFFF.DBA/FFFFA.A.B/FEE;DBA..AD;B/F/EA9B/?FFFEFFFA;FFBFFFFFA..B>.FFFFFFFFFGE.EGC-E9??CCDBA.FEFCGFGEE.EEBGBHGHHGCAGFGFEC@FEGGAFFGGGGCDFHHHHFFAF0GGDGFHHHHHF??/F//HFFFF3BEHFH@GGFB1EEEE?E3HGGGGGEE10G1GEEA2@A?BD3AAAA\tNM:i:4\tMD:Z:93A0C4A32T117\tAS:i:230\tXS:i:0\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:16398:29280\t99\tFO082262.1\t100253\t60\t237M13S\t=\t100253\t236\tGCAATCGTCGCCATGGAGAAGTATTGCGAAGAGGTTGAAAACTTCAACCAAAACTGGTGGCAGTTGATGACTGAAAACTGCACGTACGAAGAGAAGGCGAAGATGTTATCGATGGCGGACGAATCCTCGGAGGAGTATTATTTACAGGTTGGGATGCGGTTGCTCGTACCTCCGACGGGGACGGGGACGACGACGTCGTCGTCGTCTTCGGAGAGGGAGGGAATACGGACGTTCGGAGATCGGAAGAGCA\tBBBBBFABBBDAEGGGFGBG4BBGFHFEGE?FHHHHHHHHFFHHGHGFAEHHHHHHHFGEH2EHFBHDGHGHHGFHHHFHHFHFHHHGGHGFH1FG1GEGGGGHFHHFHFHHFF/GEGGGGGGHGGHHGGA@-EGCC;B;B?;BBB>BB=AB@EB=9=BFF?FEDA--DD?9.-B.;FFBEDDB=.:;ABF==FFEAFFF.ABB\tNM:i:1\tMD:Z:76C160\tAS:i:232\tXS:i:0\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:16398:29280\t147\tFO082262.1\t100253\t60\t14S236M\t=\t100253\t-236\tCGCTCTTCCGATCTGCAATCGTCGCCATGGAGAAGTATTGCGAAGAGGTTGAAAACTTCAACCAAAACTGGTGGCAGTTGATGACTGAAAACTGCACGTACGAAGAGAAGGCGAAGATGTTATCGATGGCGGACGAATCCTCGGAGGAGTATTATTTACAGGTTGGGATGCGGTTGCTCGTACCTCCGACGGGGACGGGGACGACGACGTCGTCGTCGTCTTCGGAGAGGGAGGGAATACGGACGTTCGG\t?FFFFEFFFFFFFFFFEE;BBBBFFFFFBFFEFFFFFFFFFFFFFBFFFFFFFFBFFFFEE-FFFFFFFFFFFFFFFFFGGFEGFGGGGGFGGGGGGGGFEGGGFGGGGGGHFCC0HHHGG.CDFFGGGCFGFEHHFGCGGGFHHHGHHHGDHHGHGHHHGGHGHHGGGCCGFHD@EGGHGECEE?GGGGGGGGGGHFGFGGFGGE?GGGGGGCGGGHEAGGGHHHGHGE4DFFGGEF>2AAA@BABBAA\tNM:i:1\tMD:Z:76C159\tAS:i:231\tXS:i:0\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:16195:29320\t99\tFO082266.1\t424313\t60\t250M\t=\t424557\t493\tTGTATCACATATATGAATACACAAACAAACAAACAAACGTGTGCTGAAAAGTGCAGATTTTTTCCCGCCAAAATTGTTCTCCATTGTAATGGAGCACAAAAGAGTTGTGTTTATTAAATGAATAACAACAACACCAAAAACAAAAGGTTCGCTCGGAAACCGAGACGAACGAGTGCTGCTGCGACTGCGACTGCTGCGCCTGTTTCCTCTGAAGGTGGTGTCGTTCTCGAAGAAGAGAAGGAGGAGAAGA\tBABBBBD5DBFBFGFBGGGGGGHGGHHHAGHDHHHFFHEFFHGGHHFHHHHFHHHBAFGHHHGGHHG?EEEHHHHHHHHHHHHGHGGGHHHHHCHHHHGHGBFH4FHFGHHHHHHHHHHGHHHHHHHHHGHGFF?HHHHFGHHHHHFHFFDFDFDCCGGHHGFCFGGGGGGGDGAFFCGGHHGCGGGGDGGGGGGFFGAGGEGGFGFFFGGGGGFBFFFEF/.;AAEFFE.:D;BFE.;BFF?DDFFFFF\tNM:i:0\tMD:Z:250\tAS:i:250\tXS:i:0\tRG:Z:1613\n", "M02713:279:GW200319000:1:1110:16195:29320\t147\tFO082266.1\t424557\t60\t249M\t=\t424313\t-493\tAGAAGAATGATGTTGACGATGATGATGAAATCATCATCATCGAAGAATCTTTAAAGAATAATAATAAGAAGAAGCATTCAAAGACACCGTCAAAAACGACAACGACAACATCAACGAAAAGGAAGAAGAACGATACTACCGATAATAAGAAAAGCAGTGATAAAAATGATCCGTATTCGTTCGAAAACGCGAAAGAGTTGTTCGGAAAAGTGGACGAGGAAATTTTGGTCGAGGATGATGATTATGAAG\tFBFFCBFFGBBGGGGGGGGGGGFGFGGGGGGFFGFFFFGFGHHHFHHHHHFHHEHHG0HHGFHHHHHFHBFHGBDEHEE2FGGGGGGHFHHFGAHGGGGGFFHHHCHGGGHFFGHFFHHGCEEGCGGGGF3GGFFFFDB3FA@AAA\tNM:i:2\tMD:Z:19C222G6\tAS:i:239\tXS:i:0\tRG:Z:1613\n" ] } ], "source": [ "tail 1613.sam" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Sort sam file and convert it into a bam file (compressed format) - `samtools sort`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[bam_sort_core] merging from 0 files and 8 in-memory blocks...\n" ] } ], "source": [ "samtools sort -@ 8 -o 1613.sorted.bam 1613.sam" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Check that the bam file have been created \n", "\n", "* Have a look at the filesize of the sam and bam files.\n", "* Remove the sam file " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 1.1G\n", "-rw-r--r-- 1 jovyan users 902M Oct 3 14:01 1613.sam\n", "-rw-r--r-- 1 jovyan users 167M Oct 3 14:03 1613.sorted.bam\n" ] } ], "source": [ "ls -lh\n", "rm 1613.sam" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate some stats about mapping with `samtools flagstat`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1492864 + 0 in total (QC-passed reads + QC-failed reads)\n", "1471606 + 0 primary\n", "21258 + 0 secondary\n", "0 + 0 supplementary\n", "0 + 0 duplicates\n", "0 + 0 primary duplicates\n", "1439155 + 0 mapped (96.40% : N/A)\n", "1417897 + 0 primary mapped (96.35% : N/A)\n", "1471606 + 0 paired in sequencing\n", "735803 + 0 read1\n", "735803 + 0 read2\n", "1392750 + 0 properly paired (94.64% : N/A)\n", "1413380 + 0 with itself and mate mapped\n", "4517 + 0 singletons (0.31% : N/A)\n", "16728 + 0 with mate mapped to a different chr\n", "14781 + 0 with mate mapped to a different chr (mapQ>=5)\n" ] } ], "source": [ "samtools flagstat 1613.sorted.bam" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate some stats about mapping with `samtools flagstat` and save results in the file `1613.flagstat`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "samtools flagstat 1613.sorted.bam >1613.flagstat " ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 170276\n", "-rw-r--r-- 1 jovyan users 174358429 Oct 3 14:03 1613.sorted.bam\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n" ] } ], "source": [ "ls -lrt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Display the content of the flagstat file" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1492864 + 0 in total (QC-passed reads + QC-failed reads)\n", "1471606 + 0 primary\n", "21258 + 0 secondary\n", "0 + 0 supplementary\n", "0 + 0 duplicates\n", "0 + 0 primary duplicates\n", "1439155 + 0 mapped (96.40% : N/A)\n", "1417897 + 0 primary mapped (96.35% : N/A)\n", "1471606 + 0 paired in sequencing\n", "735803 + 0 read1\n", "735803 + 0 read2\n", "1392750 + 0 properly paired (94.64% : N/A)\n", "1413380 + 0 with itself and mate mapped\n", "4517 + 0 singletons (0.31% : N/A)\n", "16728 + 0 with mate mapped to a different chr\n", "14781 + 0 with mate mapped to a different chr (mapQ>=5)\n" ] } ], "source": [ "cat 1613.flagstat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate a bam file that contains only the reads correctly paired mapped `samtools view`\n", "\n", "https://broadinstitute.github.io/picard/explain-flags.html" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "samtools view -bh -@4 -f 0x02 -o 1613.onlymappedpaired.bam 1613.sorted.bam" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 321M\n", "-rw-r--r-- 1 jovyan users 167M Oct 3 14:03 1613.sorted.bam\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n", "-rw-r--r-- 1 jovyan users 154M Oct 3 14:07 1613.onlymappedpaired.bam\n" ] } ], "source": [ "ls -lhrt" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Indexing the last bam file generated that will be used to detect SNPs" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 321M\n", "-rw-r--r-- 1 jovyan users 167M Oct 3 14:03 1613.sorted.bam\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n", "-rw-r--r-- 1 jovyan users 154M Oct 3 14:07 1613.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 46K Oct 3 14:09 1613.onlymappedpaired.bam.bai\n" ] } ], "source": [ "samtools index 1613.onlymappedpaired.bam\n", "ls -lhrt" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Let's mapping an other sample (1868) to reduce steps and tempory bam files using the pipe `|`\n", "\n", "```The pipe is a form of redirection in Linux used to connect the STDOUT of one command into the STDIN of a second command.```\n", "\n", "#### 1- Generate the bam file 1868.sorted.bam\n", "Use a single command line to perform :\n", "* mapping with `bwa-mem2 mem`\n", "* sort the mapping results with `samtools sort` \n", "* and only the file `1868.sorted.bam` is generated" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-----------------------------\n", "Executing in AVX2 mode!!\n", "-----------------------------\n", "Ref file: ../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "Entering FMI_search\n", "reference seq len = 30148641\n", "count\n", "0,\t1\n", "1,\t7834026\n", "2,\t15074321\n", "3,\t22314616\n", "4,\t30148641\n", "\n", "Reading other elements of the index from files ../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "prefix: ../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "[M::bwa_idx_load_ele] read 0 ALT contigs\n", "Done reading Index!!\n", "Reading reference genome..\n", "Binary seq file = ../REF/GCA_002220235.1_ASM222023v1_genomic.fna.0123\n", "Reference genome size: 30148640 bp\n", "Done readng reference genome !!\n", "\n", "[0000] 1: Calling process()\n", "\n", "Threads used (compute): 8\n", "Info: projected #read in a task: 529811\n", "------------------------------------------\n", "Memory pre-allocation for chaining: 1114.7223 MB\n", "Memory pre-allocation for BSW: 1916.9362 MB\n", "Memory pre-allocation for BWT: 618.5134 MB\n", "------------------------------------------\n", "No. of pipeline threads: 2\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000378, nseq: 320474\n", "\t[0000][ M::kt_pipeline] read 320474 sequences (80000378 bp)...\n", "[0000] 2. Calling mem_process_seqs.., task: 0\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000317, nseq: 320476\n", "\t[0000][ M::kt_pipeline] read 320476 sequences (80000317 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320474\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (1, 140082, 4, 2)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (215, 256, 308)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (29, 494)\n", "[M::mem_pestat] mean and std.dev: (264.14, 71.38)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 587)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320474 reads in 38.813 CPU sec, 4.866 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 1\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000454, nseq: 320474\n", "\t[0000][ M::kt_pipeline] read 320474 sequences (80000454 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320476\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (6, 139952, 4, 2)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (215, 256, 308)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (29, 494)\n", "[M::mem_pestat] mean and std.dev: (264.15, 71.74)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 587)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320476 reads in 38.959 CPU sec, 4.800 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 2\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 80000276, nseq: 320476\n", "\t[0000][ M::kt_pipeline] read 320476 sequences (80000276 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320474\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 139828, 5, 6)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (214, 255, 308)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (26, 496)\n", "[M::mem_pestat] mean and std.dev: (263.87, 72.06)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 590)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320474 reads in 38.401 CPU sec, 4.718 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 3\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 24368032, nseq: 97618\n", "\t[0000][ M::kt_pipeline] read 97618 sequences (24368032 bp)...\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 320476\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 139969, 5, 2)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (215, 256, 307)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (31, 491)\n", "[M::mem_pestat] mean and std.dev: (263.51, 71.18)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 583)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 320476 reads in 37.521 CPU sec, 4.724 real sec\n", "[0000] 2. Calling mem_process_seqs.., task: 4\n", "[0000] 3. Calling kt_for - worker_bwt\n", "[0000] read_chunk: 80000000, work_chunk_size: 0, nseq: 0\n", "[0000] 3. Calling kt_for - worker_aln\n", "Inferring insert size distribution from data, l_pac: 15074320, n: 97618\n", "[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 42657, 1, 0)\n", "[M::mem_pestat] skip orientation FF as there are not enough pairs\n", "[M::mem_pestat] analyzing insert size distribution for orientation FR...\n", "[M::mem_pestat] (25, 50, 75) percentile: (214, 256, 307)\n", "[M::mem_pestat] low and high boundaries for computing mean and std.dev: (28, 493)\n", "[M::mem_pestat] mean and std.dev: (263.48, 71.13)\n", "[M::mem_pestat] low and high boundaries for proper pairs: (1, 586)\n", "[M::mem_pestat] skip orientation RF as there are not enough pairs\n", "[M::mem_pestat] skip orientation RR as there are not enough pairs\n", "[0000] 10. Calling kt_for - worker_sam\n", "\t[0000][ M::mem_process_seqs] Processed 97618 reads in 11.981 CPU sec, 1.499 real sec\n", "[0000] read_chunk: 80000000, work_chunk_size: 0, nseq: 0\n", "[0000] Computation ends..\n", "No. of OMP threads: 8\n", "Processor is runnig @2295.055392 MHz\n", "Runtime profile:\n", "\n", "\t Time taken for main_mem function: 21.98 Sec\n", "\n", "\tIO times (sec) :\n", "\tReading IO time (reads) avg: 4.64, (4.64, 4.64)\n", "\tWriting IO time (SAM) avg: 1.86, (1.86, 1.86)\n", "\tReading IO time (Reference Genome) avg: 0.02, (0.02, 0.02)\n", "\tIndex read time avg: 0.16, (0.16, 0.16)\n", "\n", "\tOverall time (sec) (Excluding Index reading time):\n", "\tPROCESS() (Total compute time + (read + SAM) IO time) : 21.80\n", "\tMEM_PROCESS_SEQ() (Total compute time (Kernel + SAM)), avg: 20.60, (20.60, 20.60)\n", "\n", "\t SAM Processing time (sec):\n", "\t--WORKER_SAM avg: 2.51, (2.51, 2.51)\n", "\n", "\tKernels' compute time (sec):\n", "\tTotal kernel (smem+sal+bsw) time avg: 17.86, (17.86, 17.86)\n", "\t\tSMEM compute avg: 7.31, (7.37, 7.26)\n", "\t\tSAL compute avg: 0.73, (0.76, 0.69)\n", "\t\tBSW time, avg: 9.39, (9.45, 9.28)\n", "\n", "\tTotal allocs: 1640189 = 1640189 out total requests: 1975792, Rate: 0.83\n", "\n", "Important parameter settings: \n", "\tBATCH_SIZE: 512\n", "\tMAX_SEQ_LEN_REF: 256\n", "\tMAX_SEQ_LEN_QER: 128\n", "\tMAX_SEQ_LEN8: 128\n", "\tSEEDS_PER_READ: 500\n", "\tSIMD_WIDTH8 X: 32\n", "\tSIMD_WIDTH16 X: 16\n", "\tAVG_SEEDS_PER_READ: 64\n", "[bam_sort_core] merging from 0 files and 8 in-memory blocks...\n" ] } ], "source": [ "bwa-mem2 mem -R \"@RG\\tID:1868\\tSM:1868\" -t 8 ../REF/GCA_002220235.1_ASM222023v1_genomic.fna ../SHORT_READS/1868_R1.fastq.gz ../SHORT_READS/1868_R2.fastq.gz | samtools sort -@ 8 -o 1868.sorted.bam" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 489304\n", "-rw-r--r-- 1 jovyan users 174358429 Oct 3 14:03 1613.sorted.bam\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n", "-rw-r--r-- 1 jovyan users 161308863 Oct 3 14:07 1613.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 46104 Oct 3 14:09 1613.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 165320351 Oct 3 14:28 1615.sorted.bam\n" ] } ], "source": [ "ls -lrt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Filtering with samtools view and create the index" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "samtools view -bh -@8 -f 0x02 -o 1868.onlymappedpaired.bam 1868.sorted.bam\n", "samtools index 1868.onlymappedpaired.bam" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Getting some stats on mapping step - `samtools flagstat`\n", "\n", "We want the stat on the two bam files generated :\n", "* 1868.sorted.bam\n", "* 1868.onlymappedpaired.bam" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "samtools flagstat 1868.sorted.bam >1868.flagstat \n", "samtools flagstat 1868.onlymappedpaired.bam > 1868.onlymappedpaired.flagstat " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 660M\n", "-rw-r--r-- 1 jovyan users 512 Oct 3 15:42 1868.onlymappedpaired.flagstat\n", "-rw-r--r-- 1 jovyan users 519 Oct 3 15:42 1868.flagstat\n", "-rw-r--r-- 1 jovyan users 46K Oct 3 15:41 1868.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 137M Oct 3 15:41 1868.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 158M Oct 3 15:41 1868.sorted.bam\n", "-rw-r--r-- 1 jovyan users 36M Oct 3 15:34 1616_1868.vcf\n", "-rw-r--r-- 1 jovyan users 8.8M Oct 3 15:01 1616_1868.bcf\n", "-rw-r--r-- 1 jovyan users 46K Oct 3 14:09 1613.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 154M Oct 3 14:07 1613.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n", "-rw-r--r-- 1 jovyan users 167M Oct 3 14:03 1613.sorted.bam\n" ] } ], "source": [ "ls -lht" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compare the two flagstat files. What do you think about the results ?" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1393482 + 0 in total (QC-passed reads + QC-failed reads)\n", "1379518 + 0 primary\n", "0 + 0 secondary\n", "13964 + 0 supplementary\n", "0 + 0 duplicates\n", "0 + 0 primary duplicates\n", "1268617 + 0 mapped (91.04% : N/A)\n", "1254653 + 0 primary mapped (90.95% : N/A)\n", "1379518 + 0 paired in sequencing\n", "689759 + 0 read1\n", "689759 + 0 read2\n", "1241090 + 0 properly paired (89.97% : N/A)\n", "1252774 + 0 with itself and mate mapped\n", "1879 + 0 singletons (0.14% : N/A)\n", "8814 + 0 with mate mapped to a different chr\n", "7818 + 0 with mate mapped to a different chr (mapQ>=5)\n", "1249362 + 0 in total (QC-passed reads + QC-failed reads)\n", "1241090 + 0 primary\n", "0 + 0 secondary\n", "8272 + 0 supplementary\n", "0 + 0 duplicates\n", "0 + 0 primary duplicates\n", "1249362 + 0 mapped (100.00% : N/A)\n", "1241090 + 0 primary mapped (100.00% : N/A)\n", "1241090 + 0 paired in sequencing\n", "620545 + 0 read1\n", "620545 + 0 read2\n", "1241090 + 0 properly paired (100.00% : N/A)\n", "1241090 + 0 with itself and mate mapped\n", "0 + 0 singletons (0.00% : N/A)\n", "0 + 0 with mate mapped to a different chr\n", "0 + 0 with mate mapped to a different chr (mapQ>=5)\n" ] } ], "source": [ "cat 1868*.flagstat" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## SNPs calling using data from 1613 and 1868 samples - `bcftools mpileup ... | bcftools call ...`\n", "\n", "#### Print the help of `bcftools mpileup`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#PUT COMMAND" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bcftools mpileup" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "#### Print the help of `bcftools call`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#PUT COMMAND" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bcftools call" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#PUT COMMAND" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid\n", "[mpileup] 2 samples in 2 input files\n", "[mpileup] maximum number of reads per input file set to -d 250\n" ] } ], "source": [ "bcftools mpileup --threads 8 -f ../REF/GCA_002220235.1_ASM222023v1_genomic.fna *paired.bam | bcftools call --threads 8 -cv -Oz -o 1613_1868.bcf" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 638540\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n", "-rw-r--r-- 1 jovyan users 161308863 Oct 3 14:07 1613.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 46104 Oct 3 14:09 1613.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 174358429 Oct 3 14:03 1613.sorted.bam\n", "-rw-r--r-- 1 jovyan users 519 Oct 3 14:32 1615.flagstat\n", "-rw-r--r-- 1 jovyan users 143633305 Oct 3 14:30 1615.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 46184 Oct 3 14:30 1615.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 512 Oct 3 14:32 1615.onlymappedpaired.flagstat\n", "-rw-r--r-- 1 jovyan users 165320351 Oct 3 14:28 1615.sorted.bam\n", "-rw-r--r-- 1 jovyan users 9123627 Oct 3 15:01 1616_1868.bcf\n" ] } ], "source": [ "ls -l " ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [], "source": [ "bcftools convert -Ov -o 1613_1868.vcf 1613_1868.bcf\n" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total 660M\n", "-rw-r--r-- 1 jovyan users 521 Oct 3 14:06 1613.flagstat\n", "-rw-r--r-- 1 jovyan users 154M Oct 3 14:07 1613.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 46K Oct 3 14:09 1613.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 167M Oct 3 14:03 1613.sorted.bam\n", "-rw-r--r-- 1 jovyan users 519 Oct 3 14:32 1615.flagstat\n", "-rw-r--r-- 1 jovyan users 137M Oct 3 14:30 1615.onlymappedpaired.bam\n", "-rw-r--r-- 1 jovyan users 46K Oct 3 14:30 1615.onlymappedpaired.bam.bai\n", "-rw-r--r-- 1 jovyan users 512 Oct 3 14:32 1615.onlymappedpaired.flagstat\n", "-rw-r--r-- 1 jovyan users 158M Oct 3 14:28 1615.sorted.bam\n", "-rw-r--r-- 1 jovyan users 8.8M Oct 3 15:01 1616_1868.bcf\n", "-rw-r--r-- 1 jovyan users 36M Oct 3 15:34 1616_1868.vcf\n" ] } ], "source": [ "ls -lh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Have a Look to the vcf created - `head, tail`" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "##fileformat=VCFv4.2\n", "##FILTER=\n", "##bcftoolsVersion=1.13+htslib-1.13+ds\n", "##bcftoolsCommand=mpileup --threads 8 -f ../REF/GCA_002220235.1_ASM222023v1_genomic.fna 1613.onlymappedpaired.bam 1615.onlymappedpaired.bam\n", "##reference=file://../REF/GCA_002220235.1_ASM222023v1_genomic.fna\n", "##contig=\n", "##contig=\n", "##contig=\n", "##contig=\n", "##contig=\n" ] } ], "source": [ "head *vcf" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FO082259.2\t71347\t.\tG\tA\t214.22\t.\tDP=144;VDB=0.492924;SGB=29.6354;RPBZ=1.10961;BQBZ=1.78211;SCBZ=1.59266;FS=0;MQ0F=0;AF1=0.5;AC1=2;DP4=21,19,22,23;MQ=60;FQ=217.229;PV4=0.829011,1,1,1\tGT:PL\t0/0:0,120,255\t1/1:255,135,0\n", "FO082259.2\t71427\t.\tT\tG\t999\t.\tDP=157;VDB=0.619426;SGB=-1.38629;FS=0;MQ0F=0;AF1=1;AC1=4;DP4=0,0,43,47;MQ=60;FQ=-160.489\tGT:PL\t1/1:255,129,0\t1/1:255,141,0\n", "FO082259.2\t71435\t.\tG\tT\t214.22\t.\tDP=153;VDB=0.756384;SGB=29.6354;RPBZ=1.95339;BQBZ=0.568112;SCBZ=0.887406;FS=0;MQ0F=0;AF1=0.5;AC1=2;DP4=24,20,19,26;MQ=60;FQ=217.229;PV4=0.291708,0.4528,1,0.14114\tGT:PL\t0/0:0,132,255\t1/1:255,135,0\n", "FO082259.2\t71520\t.\tG\tA\t999\t.\tDP=128;VDB=0.651316;SGB=-1.38598;FS=0;MQ0F=0;AF1=1;AC1=4;DP4=0,0,39,34;MQ=60;FQ=-124.752\tGT:PL\t1/1:255,126,0\t1/1:255,93,0\n", "FO082259.2\t71540\t.\tT\tC\t999\t.\tDP=119;VDB=0.403379;SGB=24.0903;RPBZ=-1.18859;BQBZ=-0.483165;SCBZ=-1.23039;FS=0;MQ0F=0;AF1=0.5;AC1=2;DP4=16,15,17,20;MQ=60;FQ=217.229;PV4=0.807839,0.307267,1,1\tGT:PL\t1/1:255,111,0\t0/0:0,93,255\n", "FO082259.2\t71847\t.\tT\tG\t999\t.\tDP=73;VDB=0.937073;SGB=18.5451;RPBZ=-0.728589;BQBZ=-1.29364;SCBZ=-2.47615;FS=0;MQ0F=0;AF1=0.500004;AC1=2;DP4=4,13,19,10;MQ=60;FQ=217.229;PV4=0.0133826,0.0457579,1,1\tGT:PL\t1/1:255,87,0\t0/0:0,51,255\n", "FO082259.2\t71975\t.\tC\tG\t999\t.\tDP=79;VDB=0.315335;SGB=18.5451;RPBZ=-0.103735;BQBZ=-0.534962;SCBZ=-0.971325;FS=0;MQ0F=0;AF1=0.500032;AC1=2;DP4=4,10,10,19;MQ=60;FQ=217.23;PV4=1,0.249632,1,0.172793\tGT:PL\t1/1:255,87,0\t0/0:0,42,255\n", "FO082259.2\t71996\t.\tT\tG\t999\t.\tDP=88;VDB=0.965053;SGB=-1.32836;FS=0;MQ0F=0;AF1=1;AC1=4;DP4=0,0,18,32;MQ=60;FQ=-82.7539\tGT:PL\t1/1:255,99,0\t1/1:255,51,0\n", "FO082259.2\t72324\t.\tC\tT\t214.22\t.\tDP=85;VDB=0.944427;SGB=15.7726;RPBZ=0.476356;BQBZ=1.90364;SCBZ=2.99153;FS=0;MQ0F=0;AF1=0.5;AC1=2;DP4=17,10,17,8;MQ=60;FQ=217.229;PV4=0.775612,1,1,0.49283\tGT:PL\t0/0:0,81,255\t1/1:255,75,0\n", "FO082259.2\t72417\t.\tA\tC\t999\t.\tDP=90;VDB=0.69783;SGB=21.3177;RPBZ=0.133118;BQBZ=-2.21787;SCBZ=-2.50838;FS=0;MQ0F=0;AF1=0.5;AC1=2;DP4=12,9,18,15;MQ=60;FQ=217.229;PV4=1,0.0100207,1,1\tGT:PL\t1/1:255,99,0\t0/0:0,63,255\n" ] } ], "source": [ "tail 1613_1868.vcf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Download additional notebooks for exercices about mappin and SNP calling" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Go into the correct directory to download notebook - cd command\n", "cd /home/jovyan/NOTEBOOK/DAY3-4 \n", "wget https://itrop.ird.fr/trainingts/notebook/DAY3-4/Day4_Mapping_Exercice1.ipynb --no-check-certificate\n", "wget https://itrop.ird.fr/trainingts/notebook/DAY3-4/Day4_Mapping_Exercice2.ipynb --no-check-certificate" ] } ], "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": 4 }