{
"cells": [
{
"cell_type": "markdown",
"id": "6be1a832-ced2-4cfa-94bf-f23eb6558982",
"metadata": {},
"source": [
"# Formation South Green 2022 - Structural Variants Detection by using short and long reads \n",
"\n",
"# __DAY 2 : SNP calling__\n",
"\n",
"Created by C. Tranchant (DIADE-IRD), J. Orjuela (DIADE-IRD), F. Sabot (DIADE-IRD) and A. Dereeper (PHIM-IRD)\n",
"\n",
"\n",
"***\n",
"\n",
"# Table of contents\n",
"\n",
"[I - Preparing SNP calling](#data)\n",
"\n",
" * [Indexing reference with `samtools` and `GATK`](#indexref)\n",
" \n",
"[II - STEPS ANALYSYS BY SAMPLES](#gatk) \n",
"\n",
" * [Mark duplicates - gatk MarkDuplicate ](#mk)\n",
" * [Starting SNP Calling with gatk HaplotypeCaller](#HC)\n",
" * [Error produced here is normal ! RG is important in mapping step !](#BUG)đ \n",
" * [Run all the commands on each clone - `for` loop ¶](#loop)\n",
"\n",
"\n",
"[III - COMMON STEPS ANALYSES](#common)\n",
" * [Combine all the GVCF in on gvcf `gatk CombineGVCF`](#samtoolsview)\n",
" * [Compute the Genotypes - `gatk GenotypeGVCFs`](#combine)\n",
"\n",
" \n",
"\n",
"***"
]
},
{
"cell_type": "markdown",
"id": "895bacef-5993-46e5-ba43-86797d5a1ea0",
"metadata": {},
"source": [
"## __I - Preparing SNP calling__ "
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "fceac7fa-6c14-41f3-bccc-492edca10f06",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BAM_ILL MAPPING-ONT training_SV_teaching\n",
"day4-WGA-SV myFirstJupyterBook.ipynb VCF\n",
"MAPPING-ILL SV_DATA VCF_CLONES\n"
]
}
],
"source": [
"# go to work directory\n",
"cd /home/jovyan/work/\n",
"ls"
]
},
{
"cell_type": "markdown",
"id": "13a96531-2713-439a-b9dd-545e3143ad15",
"metadata": {},
"source": [
"### Go into the directory that contains all the results of the mapping performed previously with the illumina data \n",
"\n",
"Then list the content of the directory"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "fefe4d19-a627-4b8d-88bb-444918f0a0b9",
"metadata": {},
"outputs": [],
"source": [
"cd /home/jovyan/work/MAPPING-ILL"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0b3ddde3-07a2-44d4-8a06-55ba6ec2e670",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BAM dirClone11 dirClone15 dirClone19 dirClone4 dirClone8\n",
"dirClone1 dirClone12 dirClone16 dirClone2 dirClone5 dirClone9\n",
"dirClone{10} dirClone13 dirClone17 dirClone20 dirClone6 FLAGSTAT\n",
"dirClone10 dirClone14 dirClone18 dirClone3 dirClone7\n"
]
}
],
"source": [
"ls /home/jovyan/work/MAPPING-ILL"
]
},
{
"cell_type": "markdown",
"id": "85dd67e1-6696-4c37-af1e-ca381a34fc68",
"metadata": {
"tags": []
},
"source": [
"### Choose one clone to perform the SNP calling \n",
"\n",
"* We will test (almost) all the commands on a clone at first\n",
"* Let's go to declare our bash variables that will contain the path of data used in the following steps"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "a3f2f15f-a724-4e68-9eff-29d6a9995648",
"metadata": {},
"outputs": [],
"source": [
"## declare variables\n",
"i=10\n",
"REF_DIR=\"/home/jovyan/work/SV_DATA/REF/\"\n",
"REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\"\n",
"BAM=\"/home/jovyan/work/MAPPING-ILL/dirClone$i/Clone$i.SORTED.bam\""
]
},
{
"cell_type": "markdown",
"id": "8836759d-5737-4dba-ae5f-cf81a1485a8b",
"metadata": {},
"source": [
"### Indexing reference with `samtools` and `GATK` \n",
"\n",
"We will generate new index files for the reference (mandatory for the following commands)"
]
},
{
"cell_type": "markdown",
"id": "8fd685bc-62c3-44ee-aa9c-eeb05248438f",
"metadata": {},
"source": [
"#### Index the bam file just created and check that the index have been correctly created"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9a944e78-53d4-44ea-93b8-9fc19b99219f",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using GATK jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar\n",
"Running:\n",
" java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar CreateSequenceDictionary -R /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
"INFO\t2022-06-21 22:40:48\tCreateSequenceDictionary\tOutput dictionary will be written in /home/jovyan/work/SV_DATA/REF/reference.dict\n",
"22:40:48.687 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_compression.so\n",
"[Tue Jun 21 22:40:48 UTC 2022] CreateSequenceDictionary --REFERENCE /home/jovyan/work/SV_DATA/REF/reference.fasta --TRUNCATE_NAMES_AT_WHITESPACE true --NUM_SEQUENCES 2147483647 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false\n",
"[Tue Jun 21 22:40:48 UTC 2022] picard.sam.CreateSequenceDictionary done. Elapsed time: 0.00 minutes.\n",
"Runtime.totalMemory()=1141899264\n",
"To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp\n",
"picard.PicardException: /home/jovyan/work/SV_DATA/REF/reference.dict already exists. Delete this file and try again, or specify a different output file.\n",
"\tat picard.sam.CreateSequenceDictionary.doWork(CreateSequenceDictionary.java:220)\n",
"\tat picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)\n",
"\tat org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:37)\n",
"\tat org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)\n",
"\tat org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)\n",
"\tat org.broadinstitute.hellbender.Main.main(Main.java:289)\n"
]
},
{
"ename": "",
"evalue": "3",
"output_type": "error",
"traceback": []
}
],
"source": [
"cd $REF_DIR\n",
"samtools faidx $REF\n",
"gatk CreateSequenceDictionary -R $REF"
]
},
{
"cell_type": "markdown",
"id": "0b3c4903-8e01-4acc-92c1-3b6c56630735",
"metadata": {
"tags": []
},
"source": [
"#### Check that the index have been correctly created"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "799d6353-72be-44d8-b81d-a39e79d45505",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 29168\n",
"-rw-r----- 1 jovyan users 1020013 Jun 14 05:40 reference.fasta\n",
"-rw-r--r-- 1 jovyan users 255002 Jun 20 14:53 reference.fasta.pac\n",
"-rw-r--r-- 1 jovyan users 44 Jun 20 14:53 reference.fasta.ann\n",
"-rw-r--r-- 1 jovyan users 12 Jun 20 14:53 reference.fasta.amb\n",
"-rw-r--r-- 1 jovyan users 2040002 Jun 20 14:53 reference.fasta.0123\n",
"-rw-r--r-- 1 jovyan users 14280127 Jun 20 14:53 reference.fasta.bwt.8bit.32\n",
"-rw-r--r-- 1 jovyan users 12240127 Jun 20 14:53 reference.fasta.bwt.2bit.64\n",
"-rw-r--r-- 1 jovyan users 129 Jun 20 17:17 reference.dict\n",
"-rw-r--r-- 1 jovyan users 37 Jun 21 22:40 reference.fasta.fai\n"
]
}
],
"source": [
"ls -lrt"
]
},
{
"cell_type": "markdown",
"id": "9392a387-4156-4035-9fab-02fcf528bf8b",
"metadata": {},
"source": [
"## __II - STEPS ANALYSYS BY SAMPLE__ "
]
},
{
"cell_type": "markdown",
"id": "b7861d8a-56af-4ed9-8b46-3d5583931a6b",
"metadata": {
"tags": []
},
"source": [
"### Go into the directory that contains all the bam files previously generated for the clone selected"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "1eb6344c-a5f2-4e0e-a798-6f0fbfab8483",
"metadata": {},
"outputs": [],
"source": [
"# go to SR mapping results\n",
"cd /home/jovyan/work/MAPPING-ILL/dirClone${i}"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "9b436945-2f98-4662-85a2-da1dc1d904d7",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Clone10.bam clone10.gvcf.idx Clone10.SORTED.MD.bam\n",
"Clone10.flagstat Clone10.SORTED.bam Clone10.SORTED.MD.bam.bai\n",
"clone10.gvcf Clone10.SORTED.bam.bai duplicates.10.metrics\n"
]
}
],
"source": [
"ls /home/jovyan/work/MAPPING-ILL/dirClone${i}"
]
},
{
"cell_type": "markdown",
"id": "07826ea7-341e-44ae-beb8-37aa9f6d87e8",
"metadata": {},
"source": [
"### __Mark duplicates - `gatk MarkDuplicate`__ \n",
"Duplicate reads are derived from the same original physical fragment in the DNA library. There are two types of duplicates: \n",
" * PCR duplicates \n",
" * Sequencing (various optical confusions) duplicates. \n",
" \n",
"To avoid measurement bias and reduce error, only one best copy should be kept.\n",
"\n",
"#### Run gatk MarkDuplicates"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "48bc398f-3783-4616-9c4a-faad14b818dc",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
">>>>>>>>>< MarkDuplicates in Clone10 : \n",
" gatk MarkDuplicates -I /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam -M duplicates.10.metrics -O Clone10.SORTED.MD.bam\n",
"\n",
"\n",
"Using GATK jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar\n",
"Running:\n",
" java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar MarkDuplicates -I /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam -M duplicates.10.metrics -O Clone10.SORTED.MD.bam\n",
"22:41:09.222 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_compression.so\n",
"[Tue Jun 21 22:41:09 UTC 2022] MarkDuplicates --INPUT /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam --OUTPUT Clone10.SORTED.MD.bam --METRICS_FILE duplicates.10.metrics --MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP 50000 --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 --SORTING_COLLECTION_SIZE_RATIO 0.25 --TAG_DUPLICATE_SET_MEMBERS false --REMOVE_SEQUENCING_DUPLICATES false --TAGGING_POLICY DontTag --CLEAR_DT true --DUPLEX_UMI false --ADD_PG_TAG_TO_READS true --REMOVE_DUPLICATES false --ASSUME_SORTED false --DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES --PROGRAM_RECORD_ID MarkDuplicates --PROGRAM_GROUP_NAME MarkDuplicates --READ_NAME_REGEX --OPTICAL_DUPLICATE_PIXEL_DISTANCE 100 --MAX_OPTICAL_DUPLICATE_SET_SIZE 300000 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false\n",
"INFO\t2022-06-21 22:41:09\tMarkDuplicates\tStart of doWork freeMemory: 1217508896; totalMemory: 1237319680; maxMemory: 14989918208\n",
"INFO\t2022-06-21 22:41:09\tMarkDuplicates\tReading input file and constructing read end information.\n",
"INFO\t2022-06-21 22:41:09\tMarkDuplicates\tWill retain up to 54311297 data points before spilling to disk.\n",
"WARNING: BAM index file /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam.bai is older than BAM /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam\n",
"WARNING\t2022-06-21 22:41:09\tAbstractOpticalDuplicateFinderCommandLineProgram\tA field field parsed out of a read name was expected to contain an integer and did not. Read name: Reference-Clone1072030. Cause: String 'Reference-Clone1072030' did not start with a parsable number.\n",
"INFO\t2022-06-21 22:41:12\tMarkDuplicates\tRead 218229 records. 0 pairs never matched.\n",
"INFO\t2022-06-21 22:41:13\tMarkDuplicates\tAfter buildSortedReadEndLists freeMemory: 1240047976; totalMemory: 1722810368; maxMemory: 14989918208\n",
"INFO\t2022-06-21 22:41:13\tMarkDuplicates\tWill retain up to 468434944 duplicate indices before spilling to disk.\n",
"INFO\t2022-06-21 22:41:15\tMarkDuplicates\tTraversing read pair information and detecting duplicates.\n",
"INFO\t2022-06-21 22:41:15\tMarkDuplicates\tTraversing fragment information and detecting duplicates.\n",
"INFO\t2022-06-21 22:41:15\tMarkDuplicates\tSorting list of duplicate records.\n",
"INFO\t2022-06-21 22:41:16\tMarkDuplicates\tAfter generateDuplicateIndexes freeMemory: 1697721576; totalMemory: 5470420992; maxMemory: 14989918208\n",
"INFO\t2022-06-21 22:41:16\tMarkDuplicates\tMarking 136 records as duplicates.\n",
"INFO\t2022-06-21 22:41:16\tMarkDuplicates\tFound 0 optical duplicate clusters.\n",
"WARNING: BAM index file /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam.bai is older than BAM /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam\n",
"INFO\t2022-06-21 22:41:16\tMarkDuplicates\tReads are assumed to be ordered by: coordinate\n",
"INFO\t2022-06-21 22:41:17\tMarkDuplicates\tWriting complete. Closing input iterator.\n",
"INFO\t2022-06-21 22:41:17\tMarkDuplicates\tDuplicate Index cleanup.\n",
"INFO\t2022-06-21 22:41:17\tMarkDuplicates\tGetting Memory Stats.\n",
"INFO\t2022-06-21 22:41:17\tMarkDuplicates\tBefore output close freeMemory: 5443958392; totalMemory: 5470420992; maxMemory: 14989918208\n",
"INFO\t2022-06-21 22:41:17\tMarkDuplicates\tClosed outputs. Getting more Memory Stats.\n",
"INFO\t2022-06-21 22:41:17\tMarkDuplicates\tAfter output close freeMemory: 5446297336; totalMemory: 5470420992; maxMemory: 14989918208\n",
"[Tue Jun 21 22:41:17 UTC 2022] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 0.15 minutes.\n",
"Runtime.totalMemory()=5470420992\n",
"Tool returned:\n",
"0\n"
]
}
],
"source": [
"echo -e \"\\n>>>>>>>>>< MarkDuplicates in Clone$i : \\n gatk MarkDuplicates -I $BAM -M duplicates.$i.metrics -O Clone$i.SORTED.MD.bam\\n\\n\";\n",
"gatk MarkDuplicates -I $BAM -M duplicates.$i.metrics -O Clone$i.SORTED.MD.bam;"
]
},
{
"cell_type": "markdown",
"id": "7a3148e8-1f37-4b5c-afa6-aa9120b4a39b",
"metadata": {},
"source": [
"#### Check that a bam file have been correctly created"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "db219f2c-fdfe-463f-a804-4c47abf2d70b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 79M\n",
"-rw-r--r-- 1 jovyan users 2.8K Jun 21 22:11 Clone10.SORTED.bam.bai\n",
"-rw-r--r-- 1 jovyan users 39M Jun 21 22:13 Clone10.bam\n",
"-rw-r--r-- 1 jovyan users 381 Jun 21 22:13 Clone10.flagstat\n",
"-rw-r--r-- 1 jovyan users 18M Jun 21 22:13 Clone10.SORTED.bam\n",
"-rw-r--r-- 1 jovyan users 218 Jun 21 22:14 clone10.gvcf.idx\n",
"-rw-r--r-- 1 jovyan users 13K Jun 21 22:14 clone10.gvcf\n",
"-rw-r--r-- 1 jovyan users 2.8K Jun 21 22:36 Clone10.SORTED.MD.bam.bai\n",
"-rw-r--r-- 1 jovyan users 23M Jun 21 22:41 Clone10.SORTED.MD.bam\n",
"-rw-r--r-- 1 jovyan users 3.4K Jun 21 22:41 duplicates.10.metrics\n"
]
}
],
"source": [
"ls -lrth"
]
},
{
"cell_type": "markdown",
"id": "da96d3e6-3155-40bc-810e-bdcd89eed2f3",
"metadata": {},
"source": [
"#### Index the bam file just created and check that the index have been correctly created"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "aaed6229-03c7-4946-9355-bec5eb587592",
"metadata": {},
"outputs": [],
"source": [
"samtools index Clone$i.SORTED.MD.bam;"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "6aafebd6-bc19-4494-8a57-197b68c1e8e0",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 79M\n",
"-rw-r--r-- 1 jovyan users 2.8K Jun 21 22:11 Clone10.SORTED.bam.bai\n",
"-rw-r--r-- 1 jovyan users 39M Jun 21 22:13 Clone10.bam\n",
"-rw-r--r-- 1 jovyan users 381 Jun 21 22:13 Clone10.flagstat\n",
"-rw-r--r-- 1 jovyan users 18M Jun 21 22:13 Clone10.SORTED.bam\n",
"-rw-r--r-- 1 jovyan users 218 Jun 21 22:14 clone10.gvcf.idx\n",
"-rw-r--r-- 1 jovyan users 13K Jun 21 22:14 clone10.gvcf\n",
"-rw-r--r-- 1 jovyan users 23M Jun 21 22:41 Clone10.SORTED.MD.bam\n",
"-rw-r--r-- 1 jovyan users 3.4K Jun 21 22:41 duplicates.10.metrics\n",
"-rw-r--r-- 1 jovyan users 2.8K Jun 21 22:41 Clone10.SORTED.MD.bam.bai\n"
]
}
],
"source": [
"# check MD bam\n",
"ls -lrth"
]
},
{
"cell_type": "markdown",
"id": "ddc481e5-9697-4106-a145-bb7d6183a4a0",
"metadata": {},
"source": [
"## Starting SNP Calling with `gatk HaplotypeCaller` \n",
"\n",
"HaplotypeCaller is used to call potential variant sites per sample and save results in GVCF format. With GVCF, it provides variant sites, and groups non-variant sites into blocks during the calling process based on genotype quality. \n",
" \n",
"This is a way of compressing the VCF file without losing any sites in order to do joint analysis in subsequent steps.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "1e0952cc-2a7f-4e0a-bc17-6bfc4d367e13",
"metadata": {},
"outputs": [],
"source": [
"# change BAM file for MD one\n",
"BAM=\"/home/jovyan/work/MAPPING-ILL/dirClone$i/Clone$i.SORTED.MD.bam\""
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5e8e1eb-3d05-40a5-9957-021f83ed9965",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# lauch GATK HaplotypeCaller\n",
"echo -e \"\\nCalling Clone$i \\n gatk HaplotypeCaller -I $BAM -O /home/jovyan/work/MAPPING-ILL/dirClone$i/clone${i}.gvcf -R $REF -ERC GVCF\"\n",
"gatk HaplotypeCaller -I $BAM -O /home/jovyan/work/MAPPING-ILL/dirClone$i/clone${i}.gvcf -R $REF -ERC GVCF"
]
},
{
"cell_type": "markdown",
"id": "86d95ae1-5e20-4949-bc4f-a95fb3bb8cc7",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"### Error produced here is normal WHY ?? ! đ "
]
},
{
"cell_type": "markdown",
"id": "2851f135-79da-4868-aab4-6a8a1c956cc6",
"metadata": {
"tags": []
},
"source": [
"### If you have found error, now you have to relaunch mapping with all samples ! BUT ...\n",
"\n",
"-R \"@RG\\tID:Clone${i}\\tSM:Clone${i}\"\n",
"\n",
"Don't worry, correct bam files with correct are available for downloading."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13993511-2f8a-4ebd-8b2c-dd0e82086632",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"cd /home/jovyan/work/\n",
"wget --no-check-certificat -rm -nH --cut-dirs=1 --reject=\"index.html*\" wget https://itrop.ird.fr/sv-training/BAM_ILL.tar.gz\n",
"tar zxvf BAM_ILL.tar.gz\n",
"rm BAM_ILL.tar.gz\n"
]
},
{
"cell_type": "markdown",
"id": "74813193-5d58-4a21-ac39-02233366c7f4",
"metadata": {},
"source": [
"#### Store the path of this new directory in the bash variable `BAM_ILL` "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "39ce7d41-bdbc-40f9-8fbe-09eda8ce5f85",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"BAM_ILL=\"/home/jovyan/work/BAM_ILL\""
]
},
{
"cell_type": "markdown",
"id": "28bec3c7-6a8b-4652-93b8-ac5e3b2396e0",
"metadata": {},
"source": [
"#### List the content of the directory `BAM_ILL`"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "fe2437a0-5ffa-47eb-b6b0-efa4e58ad61f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 472K\n",
"-rw-r--r-- 1 jovyan users 1.3K Jun 20 14:44 README.md\n",
"-rw-r--r-- 1 jovyan users 12K Jun 20 14:44 Day3b_Variants_structuraux_CORR.ipynb\n",
"-rw-r--r-- 1 jovyan users 19K Jun 20 17:03 Day1_Mapping_Practice_EMPTY.ipynb\n",
"-rw-r--r-- 1 jovyan users 72K Jun 20 20:21 Day2a_Mapping_analysis_CORR.ipynb\n",
"-rw-r--r-- 1 jovyan users 11K Jun 20 20:23 Day2a_Mapping_analysis-EMPTY.ipynb\n",
"-rw-r--r-- 1 jovyan users 20K Jun 20 21:35 Day2c_SNP_analysis_EMPTY.ipynb\n",
"-rw-r--r-- 1 jovyan users 144K Jun 20 21:40 Day2c_SNP_analysis_CORR.ipynb\n",
"-rw-r--r-- 1 jovyan users 15K Jun 20 21:59 Day4_SV_genome_CORR.ipynb\n",
"-rw-r--r-- 1 jovyan users 15K Jun 21 05:59 Day4_SV_genome_CORR-Copy1.ipynb\n",
"-rw-r--r-- 1 jovyan users 61K Jun 21 21:29 Day1_Mapping_Practice_CORR.ipynb\n",
"-rw-r--r-- 1 jovyan users 63K Jun 21 22:38 Day2b_Appel_variants_SNP_CORR.ipynb\n",
"-rw-r--r-- 1 jovyan users 28K Jun 21 22:52 Day2b_Appel_variants_SNP_EMPTY.ipynb\n"
]
}
],
"source": [
"ls -lrth $BAM_ILL"
]
},
{
"cell_type": "markdown",
"id": "d6649e9f-b627-4b8b-b1dc-b4ee6098a736",
"metadata": {
"tags": []
},
"source": [
"#### Create the subdirectory `VCF` into the directory `work`\n",
" \n",
"We will save all the vcf generated into this directory"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "fc1f5a57-96ff-412c-a866-f536314c1f36",
"metadata": {},
"outputs": [],
"source": [
"# create a work directory to SNP calling results\n",
"mkdir -p /home/jovyan/work/VCF\n",
"VCF_DIR=\"/home/jovyan/work/VCF\""
]
},
{
"cell_type": "markdown",
"id": "46a64c92-a390-41e2-9780-ae1739340cd1",
"metadata": {},
"source": [
"#### Run gatk markduplicates and haplotypecaller on the new bam of the clone selected \n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c0c2d255-5ea7-48b6-b3ad-0618937b149a",
"metadata": {},
"outputs": [],
"source": [
"## declare variables\n",
"i=10\n",
"REF_DIR=\"/home/jovyan/work/SV_DATA/REF/\"\n",
"REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\"\n",
"BAM=\"/home/jovyan/work/BAM_ILL/Clone$i.SORTED.bam\""
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "396746ec-919a-4129-88c9-b3073ade51fe",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
">>>>>>>>>< MarkDuplicates in Clone10 : \n",
" gatk MarkDuplicates -I /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam -M duplicates.10.metrics -O Clone10.SORTED.MD.bam\n",
"\n",
"\n",
"Using GATK jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar\n",
"Running:\n",
" java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar MarkDuplicates -I /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam -M duplicates.10.metrics -O /home/jovyan/work/BAM_ILL/Clone10.SORTED.MD.bam\n",
"23:02:12.847 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_compression.so\n",
"[Tue Jun 21 23:02:12 UTC 2022] MarkDuplicates --INPUT /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam --OUTPUT /home/jovyan/work/BAM_ILL/Clone10.SORTED.MD.bam --METRICS_FILE duplicates.10.metrics --MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP 50000 --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 --SORTING_COLLECTION_SIZE_RATIO 0.25 --TAG_DUPLICATE_SET_MEMBERS false --REMOVE_SEQUENCING_DUPLICATES false --TAGGING_POLICY DontTag --CLEAR_DT true --DUPLEX_UMI false --ADD_PG_TAG_TO_READS true --REMOVE_DUPLICATES false --ASSUME_SORTED false --DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES --PROGRAM_RECORD_ID MarkDuplicates --PROGRAM_GROUP_NAME MarkDuplicates --READ_NAME_REGEX --OPTICAL_DUPLICATE_PIXEL_DISTANCE 100 --MAX_OPTICAL_DUPLICATE_SET_SIZE 300000 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false\n",
"INFO\t2022-06-21 23:02:13\tMarkDuplicates\tStart of doWork freeMemory: 1062809128; totalMemory: 1079508992; maxMemory: 14989918208\n",
"INFO\t2022-06-21 23:02:13\tMarkDuplicates\tReading input file and constructing read end information.\n",
"INFO\t2022-06-21 23:02:13\tMarkDuplicates\tWill retain up to 54311297 data points before spilling to disk.\n",
"WARNING: BAM index file /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam.bai is older than BAM /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam\n",
"WARNING\t2022-06-21 23:02:13\tAbstractOpticalDuplicateFinderCommandLineProgram\tA field field parsed out of a read name was expected to contain an integer and did not. Read name: Reference-Clone1072030. Cause: String 'Reference-Clone1072030' did not start with a parsable number.\n",
"INFO\t2022-06-21 23:02:16\tMarkDuplicates\tRead 218229 records. 0 pairs never matched.\n",
"INFO\t2022-06-21 23:02:17\tMarkDuplicates\tAfter buildSortedReadEndLists freeMemory: 990260360; totalMemory: 1469054976; maxMemory: 14989918208\n",
"INFO\t2022-06-21 23:02:17\tMarkDuplicates\tWill retain up to 468434944 duplicate indices before spilling to disk.\n",
"INFO\t2022-06-21 23:02:18\tMarkDuplicates\tTraversing read pair information and detecting duplicates.\n",
"INFO\t2022-06-21 23:02:18\tMarkDuplicates\tTraversing fragment information and detecting duplicates.\n",
"INFO\t2022-06-21 23:02:19\tMarkDuplicates\tSorting list of duplicate records.\n",
"INFO\t2022-06-21 23:02:19\tMarkDuplicates\tAfter generateDuplicateIndexes freeMemory: 1453721656; totalMemory: 5222432768; maxMemory: 14989918208\n",
"INFO\t2022-06-21 23:02:19\tMarkDuplicates\tMarking 136 records as duplicates.\n",
"INFO\t2022-06-21 23:02:19\tMarkDuplicates\tFound 0 optical duplicate clusters.\n",
"WARNING: BAM index file /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam.bai is older than BAM /home/jovyan/work/MAPPING-ILL/dirClone10/Clone10.SORTED.bam\n",
"INFO\t2022-06-21 23:02:19\tMarkDuplicates\tReads are assumed to be ordered by: coordinate\n",
"INFO\t2022-06-21 23:02:20\tMarkDuplicates\tWriting complete. Closing input iterator.\n",
"INFO\t2022-06-21 23:02:20\tMarkDuplicates\tDuplicate Index cleanup.\n",
"INFO\t2022-06-21 23:02:20\tMarkDuplicates\tGetting Memory Stats.\n",
"INFO\t2022-06-21 23:02:20\tMarkDuplicates\tBefore output close freeMemory: 1463182176; totalMemory: 5223481344; maxMemory: 14989918208\n",
"INFO\t2022-06-21 23:02:21\tMarkDuplicates\tClosed outputs. Getting more Memory Stats.\n",
"INFO\t2022-06-21 23:02:21\tMarkDuplicates\tAfter output close freeMemory: 5790348560; totalMemory: 5821693952; maxMemory: 14989918208\n",
"[Tue Jun 21 23:02:21 UTC 2022] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 0.14 minutes.\n",
"Runtime.totalMemory()=5821693952\n",
"Tool returned:\n",
"0\n",
"Using GATK jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar\n",
"Running:\n",
" java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xmx20g -jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar HaplotypeCaller -ERC GVCF -I /home/jovyan/work/BAM_ILL/Clone10.SORTED.MD.bam -O /home/jovyan/work/VCF/clone10.gvcf -R /home/jovyan/work/SV_DATA/REF/reference.fasta\n",
"23:02:25.749 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_compression.so\n",
"23:02:25.938 INFO HaplotypeCaller - ------------------------------------------------------------\n",
"23:02:25.939 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.2.6.1\n",
"23:02:25.939 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/\n",
"23:02:25.948 INFO HaplotypeCaller - Initializing engine\n",
"23:02:26.380 INFO HaplotypeCaller - Done initializing engine\n",
"23:02:26.384 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled\n",
"23:02:26.391 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to 0.0 for reference-model confidence output\n",
"23:02:26.391 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output\n",
"23:02:26.415 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_utils.so\n",
"23:02:26.417 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so\n",
"23:02:26.438 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM\n",
"23:02:26.443 INFO IntelPairHmm - Available threads: 16\n",
"23:02:26.444 INFO IntelPairHmm - Requested threads: 4\n",
"23:02:26.444 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation\n",
"23:02:26.450 WARN GATKVariantContextUtils - Can't determine output variant file format from output file extension \"gvcf\". Defaulting to VCF.\n",
"23:02:26.501 INFO ProgressMeter - Starting traversal\n",
"23:02:26.502 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute\n",
"23:02:27.333 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position Reference:58 and possibly subsequent; at least 10 samples must have called genotypes\n",
"23:02:37.755 INFO ProgressMeter - Reference:5531 0.2 30 160.0\n",
"23:02:49.154 INFO ProgressMeter - Reference:11476 0.4 60 158.9\n",
"23:02:50.142 WARN DepthPerSampleHC - Annotation will not be calculated at position Reference:12293 and possibly subsequent; genotype for sample Clone10 is not called\n",
"23:02:50.143 WARN StrandBiasBySample - Annotation will not be calculated at position Reference:12293 and possibly subsequent; genotype for sample Clone10 is not called\n",
"23:03:00.297 INFO ProgressMeter - Reference:41847 0.6 180 319.6\n",
"23:03:10.615 INFO ProgressMeter - Reference:54347 0.7 240 326.4\n",
"23:03:22.138 INFO ProgressMeter - Reference:61384 0.9 280 302.0\n",
"23:03:32.426 INFO ProgressMeter - Reference:69383 1.1 320 291.2\n",
"23:03:42.716 INFO ProgressMeter - Reference:78459 1.3 370 291.3\n",
"23:03:54.930 INFO ProgressMeter - Reference:87547 1.5 420 285.0\n",
"23:04:06.829 INFO ProgressMeter - Reference:96058 1.7 470 281.1\n",
"23:04:18.292 INFO ProgressMeter - Reference:101304 1.9 500 268.4\n",
"23:04:31.152 INFO ProgressMeter - Reference:109970 2.1 550 264.7\n",
"23:04:42.162 INFO ProgressMeter - Reference:117189 2.3 590 260.9\n",
"23:04:52.770 INFO ProgressMeter - Reference:124698 2.4 630 258.4\n",
"23:05:03.454 INFO ProgressMeter - Reference:133648 2.6 680 260.0\n",
"23:05:15.965 INFO ProgressMeter - Reference:141299 2.8 720 254.9\n",
"23:05:25.967 INFO ProgressMeter - Reference:146913 3.0 750 250.7\n",
"23:05:36.794 INFO ProgressMeter - Reference:152244 3.2 780 245.9\n",
"23:05:47.845 INFO ProgressMeter - Reference:159190 3.4 820 244.4\n",
"23:06:00.601 INFO ProgressMeter - Reference:166631 3.6 860 241.0\n",
"23:06:10.831 INFO ProgressMeter - Reference:175186 3.7 910 243.4\n",
"23:06:23.326 INFO ProgressMeter - Reference:198991 3.9 1010 255.9\n",
"23:06:35.665 INFO ProgressMeter - Reference:208340 4.2 1060 255.3\n",
"23:06:47.719 INFO ProgressMeter - Reference:216500 4.4 1110 255.0\n",
"23:06:59.265 INFO ProgressMeter - Reference:232906 4.5 1190 261.8\n",
"23:07:10.958 INFO ProgressMeter - Reference:240444 4.7 1230 259.4\n",
"23:07:23.886 INFO ProgressMeter - Reference:252382 5.0 1300 262.3\n",
"23:07:34.380 INFO ProgressMeter - Reference:263011 5.1 1360 265.0\n",
"23:07:47.609 INFO ProgressMeter - Reference:269942 5.4 1400 261.6\n",
"23:07:59.221 INFO ProgressMeter - Reference:291661 5.5 1490 268.7\n",
"23:08:09.702 INFO ProgressMeter - Reference:301679 5.7 1550 271.0\n",
"23:08:20.339 INFO ProgressMeter - Reference:310885 5.9 1600 271.3\n",
"23:08:31.705 INFO ProgressMeter - Reference:316082 6.1 1630 267.8\n",
"23:08:43.216 INFO ProgressMeter - Reference:398244 6.3 1930 307.4\n",
"23:08:54.399 INFO ProgressMeter - Reference:410802 6.5 1990 307.8\n",
"23:09:06.314 INFO ProgressMeter - Reference:428714 6.7 2070 310.6\n",
"23:09:17.131 INFO ProgressMeter - Reference:437066 6.8 2120 309.8\n",
"23:09:30.402 INFO ProgressMeter - Reference:447863 7.1 2180 308.6\n",
"23:09:41.233 INFO ProgressMeter - Reference:457345 7.2 2230 307.8\n",
"23:09:51.426 INFO ProgressMeter - Reference:463145 7.4 2260 304.8\n",
"23:10:03.888 INFO ProgressMeter - Reference:470429 7.6 2300 301.7\n",
"23:10:14.316 INFO ProgressMeter - Reference:479231 7.8 2350 301.4\n",
"23:10:27.122 INFO ProgressMeter - Reference:501284 8.0 2440 304.6\n",
"23:10:37.930 INFO ProgressMeter - Reference:507169 8.2 2470 301.6\n",
"23:10:49.189 INFO ProgressMeter - Reference:530068 8.4 2560 305.6\n",
"23:10:59.332 INFO ProgressMeter - Reference:535643 8.5 2590 303.0\n",
"23:11:13.849 INFO ProgressMeter - Reference:542354 8.8 2630 299.2\n",
"23:11:24.199 INFO ProgressMeter - Reference:552817 9.0 2680 299.1\n",
"23:11:36.543 INFO ProgressMeter - Reference:559682 9.2 2720 296.7\n",
"23:11:47.469 INFO ProgressMeter - Reference:567556 9.3 2760 295.2\n",
"23:12:00.149 INFO ProgressMeter - Reference:586874 9.6 2840 297.0\n",
"23:12:11.477 INFO ProgressMeter - Reference:594049 9.7 2880 295.4\n",
"23:12:26.125 INFO ProgressMeter - Reference:601302 10.0 2920 292.2\n",
"23:12:39.392 INFO ProgressMeter - Reference:608767 10.2 2960 289.8\n",
"23:12:51.490 INFO ProgressMeter - Reference:613966 10.4 2990 287.0\n",
"23:13:05.230 INFO ProgressMeter - Reference:620334 10.6 3030 284.6\n",
"23:13:15.738 INFO ProgressMeter - Reference:623917 10.8 3050 281.9\n",
"23:13:29.676 INFO ProgressMeter - Reference:629293 11.1 3080 278.7\n",
"23:13:41.130 INFO ProgressMeter - Reference:634683 11.2 3110 276.6\n",
"23:13:51.990 INFO ProgressMeter - Reference:641202 11.4 3150 275.7\n",
"23:14:04.294 INFO ProgressMeter - Reference:647341 11.6 3180 273.4\n",
"23:14:15.797 INFO ProgressMeter - Reference:654590 11.8 3220 272.4\n",
"23:14:26.635 INFO ProgressMeter - Reference:660161 12.0 3250 270.8\n",
"23:14:39.336 INFO ProgressMeter - Reference:666011 12.2 3280 268.5\n",
"23:14:51.880 INFO ProgressMeter - Reference:680860 12.4 3350 269.7\n",
"23:15:05.761 INFO ProgressMeter - Reference:686407 12.7 3380 267.1\n",
"23:15:18.146 INFO ProgressMeter - Reference:691997 12.9 3410 265.1\n",
"23:15:28.157 INFO ProgressMeter - Reference:696933 13.0 3440 264.1\n",
"23:15:40.491 INFO ProgressMeter - Reference:704214 13.2 3480 263.0\n",
"23:15:52.185 INFO ProgressMeter - Reference:717990 13.4 3540 263.6\n",
"23:16:02.431 INFO ProgressMeter - Reference:728953 13.6 3600 264.7\n",
"23:16:15.307 INFO ProgressMeter - Reference:737246 13.8 3650 264.2\n",
"23:16:27.857 INFO ProgressMeter - Reference:746747 14.0 3700 263.9\n",
"23:16:41.021 INFO ProgressMeter - Reference:752794 14.2 3730 261.9\n",
"23:16:51.398 INFO ProgressMeter - Reference:763704 14.4 3790 262.9\n",
"23:17:01.490 INFO ProgressMeter - Reference:772720 14.6 3840 263.3\n",
"23:17:14.723 INFO ProgressMeter - Reference:781816 14.8 3890 262.8\n",
"23:17:27.915 INFO ProgressMeter - Reference:793881 15.0 3960 263.6\n",
"23:17:39.860 INFO ProgressMeter - Reference:802742 15.2 4010 263.4\n",
"23:17:51.732 INFO ProgressMeter - Reference:810162 15.4 4050 262.6\n",
"23:18:05.191 INFO ProgressMeter - Reference:835275 15.6 4160 265.9\n",
"23:18:15.353 INFO ProgressMeter - Reference:842464 15.8 4200 265.6\n",
"23:18:25.990 INFO ProgressMeter - Reference:847680 16.0 4230 264.5\n",
"23:18:37.396 INFO ProgressMeter - Reference:873284 16.2 4340 268.2\n",
"23:18:49.981 INFO ProgressMeter - Reference:882322 16.4 4390 267.8\n",
"23:19:03.530 INFO ProgressMeter - Reference:889700 16.6 4430 266.6\n",
"23:19:14.537 INFO ProgressMeter - Reference:896349 16.8 4470 266.1\n",
"23:19:26.190 INFO ProgressMeter - Reference:901724 17.0 4500 264.8\n",
"23:19:37.339 INFO ProgressMeter - Reference:909572 17.2 4540 264.3\n",
"23:19:47.775 INFO ProgressMeter - Reference:918168 17.4 4590 264.5\n",
"23:20:01.239 INFO ProgressMeter - Reference:925823 17.6 4630 263.4\n",
"23:20:13.048 INFO ProgressMeter - Reference:950830 17.8 4730 266.1\n",
"23:20:25.546 INFO ProgressMeter - Reference:956619 18.0 4760 264.7\n",
"23:20:35.801 INFO ProgressMeter - Reference:963096 18.2 4800 264.4\n",
"23:20:48.487 INFO ProgressMeter - Reference:973230 18.4 4860 264.6\n",
"23:21:01.771 INFO ProgressMeter - Reference:980792 18.6 4900 263.6\n",
"23:21:14.437 INFO ProgressMeter - Reference:988343 18.8 4940 262.8\n",
"23:21:19.465 WARN IntelInflater - Zero Bytes Written : 0\n",
"23:21:20.629 INFO HaplotypeCaller - 21683 read(s) filtered by: MappingQualityReadFilter \n",
"0 read(s) filtered by: MappingQualityAvailableReadFilter \n",
"0 read(s) filtered by: MappedReadFilter \n",
"132 read(s) filtered by: NotSecondaryAlignmentReadFilter \n",
"120 read(s) filtered by: NotDuplicateReadFilter \n",
"0 read(s) filtered by: PassesVendorQualityCheckReadFilter \n",
"0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter \n",
"0 read(s) filtered by: GoodCigarReadFilter \n",
"0 read(s) filtered by: WellformedReadFilter \n",
"21935 total reads filtered\n",
"23:21:20.629 INFO ProgressMeter - Reference:1017872 18.9 5067 268.1\n",
"23:21:20.630 INFO ProgressMeter - Traversal complete. Processed 5067 total regions in 18.9 minutes.\n",
"23:21:20.658 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.177231199\n",
"23:21:20.659 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 94.42924324200001\n",
"23:21:20.659 INFO SmithWatermanAligner - Total compute time in java Smith-Waterman : 924.92 sec\n",
"23:21:20.659 INFO HaplotypeCaller - Shutting down engine\n",
"[June 21, 2022 11:21:20 PM UTC] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 18.92 minutes.\n",
"Runtime.totalMemory()=1074266112\n"
]
}
],
"source": [
"echo -e \"\\n>>>>>>>>>< MarkDuplicates in Clone$i : \\n gatk MarkDuplicates -I $BAM -M duplicates.$i.metrics -O Clone$i.SORTED.MD.bam\\n\\n\";\n",
"gatk MarkDuplicates -I $BAM -M duplicates.$i.metrics -O /home/jovyan/work/BAM_ILL/Clone$i.SORTED.MD.bam\n",
"\n",
"samtools index /home/jovyan/work/BAM_ILL/Clone$i.SORTED.MD.bam\n",
"\n",
"gatk HaplotypeCaller --native-pair-hmm-threads 4 --java-options \"-Xmx20g\" -ERC GVCF -I /home/jovyan/work/BAM_ILL/Clone${i}.SORTED.MD.bam -O /home/jovyan/work/VCF/clone${i}.gvcf -R $REF"
]
},
{
"cell_type": "markdown",
"id": "2cad5bbc-ccbf-4bea-9acf-c86f78d02838",
"metadata": {},
"source": [
"#### Check that the gvcf files have been created"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "0894d2f2-fdc3-4e86-a259-e2a23b8fe062",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"clone10.gvcf Clone13.g.vcf.idx Clone18.g.vcf density.csv\n",
"Clone10.g.vcf Clone14.g.vcf Clone18.g.vcf.idx genome.txt\n",
"clone10.gvcf.idx Clone14.g.vcf.idx Clone19.g.vcf output.vcf\n",
"Clone10.g.vcf.idx Clone15.g.vcf Clone19.g.vcf.idx output.vcf.idx\n",
"Clone11.g.vcf Clone15.g.vcf.idx Clone1.g.vcf rawSNP.vcf\n",
"Clone11.g.vcf.idx Clone16.g.vcf Clone1.g.vcf.idx rawSNP.vcf.idx\n",
"Clone12.g.vcf Clone16.g.vcf.idx Clone20.g.vcf\n",
"Clone12.g.vcf.idx Clone17.g.vcf Clone20.g.vcf.idx\n",
"Clone13.g.vcf Clone17.g.vcf.idx Clone2.g.vcf\n"
]
}
],
"source": [
"ls $VCF_DIR/"
]
},
{
"cell_type": "markdown",
"id": "b75d535f-de6b-495a-9b0b-22a23bfb34f0",
"metadata": {},
"source": [
"#### Have a look at the gvfc file created"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "bc0acfeb-915e-476e-89f1-3e95ac27377a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##fileformat=VCFv4.2\n",
"##ALT=\n",
"##FILTER=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"Reference\t996185\t.\tG\tT,\t255.97\t.\tDP=6;ExcessHet=0.0000;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQandDP=21600,6\tGT:AD:DP:GQ:PGT:PID:PL:PS:SB\t1|1:0,6,0:6:18:0|1:996181_C_T:270,18,0,270,18,270:996181:0,0,0,6\n",
"Reference\t996186\t.\tT\t\t.\t.\tEND=996192\tGT:DP:GQ:MIN_DP:PL\t0/0:6:18:6:0,18,233\n",
"Reference\t996193\t.\tT\t\t.\t.\tEND=996198\tGT:DP:GQ:MIN_DP:PL\t0/0:6:12:6:0,12,180\n",
"Reference\t996199\t.\tA\t\t.\t.\tEND=996202\tGT:DP:GQ:MIN_DP:PL\t0/0:6:9:6:0,9,135\n",
"Reference\t996203\t.\tC\tT,\t255.97\t.\tDP=6;ExcessHet=0.0000;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQandDP=21600,6\tGT:AD:DP:GQ:PGT:PID:PL:PS:SB\t1|1:0,6,0:6:18:0|1:996181_C_T:270,18,0,270,18,270:996181:0,0,0,6\n",
"Reference\t996204\t.\tG\t\t.\t.\tEND=996216\tGT:DP:GQ:MIN_DP:PL\t0/0:3:9:3:0,9,114\n",
"Reference\t996217\t.\tC\t\t.\t.\tEND=996222\tGT:DP:GQ:MIN_DP:PL\t0/0:3:3:3:0,3,45\n",
"Reference\t996223\t.\tA\t\t.\t.\tEND=996227\tGT:DP:GQ:MIN_DP:PL\t0/0:3:0:3:0,0,0\n",
"Reference\t996228\t.\tG\tC,\t121.83\t.\tDP=3;ExcessHet=0.0000;MLEAC=1,0;MLEAF=0.500,0.00;RAW_MQandDP=10800,3\tGT:AD:DP:GQ:PGT:PID:PL:PS:SB\t1|1:0,3,0:3:9:0|1:996181_C_T:135,9,0,135,9,135:996181:0,0,0,3\n",
"Reference\t996229\t.\tA\t\t.\t.\tEND=1020001\tGT:DP:GQ:MIN_DP:PL\t0/0:0:0:0:0,0,0\n"
]
}
],
"source": [
"head $VCF_DIR/clone${i}.gvcf"
]
},
{
"cell_type": "markdown",
"id": "2a1a59ea-cd90-4533-a33d-17e1db77a17a",
"metadata": {},
"source": [
"## Run all the commands on each clone `for loop` "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f16f3738-431c-44c7-a4f6-e4493f92a232",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"cd $BAM_ILL\n",
"for i in {1..2}\n",
" do\n",
" samtools index Clone$i.SORTED.bam;\n",
" echo -e \"\\nMarkDuplicates in Clone$i\";\n",
" gatk MarkDuplicates -I Clone$i.SORTED.bam -M duplicates.$i.metrics -O Clone$i.SORTED.MD.bam;\n",
" samtools index Clone$i.SORTED.MD.bam;\n",
" echo -e \"\\n\\n>>>>>>>>>>>> Calling Clone$i\";\n",
" gatk --java-options \"-Xmx20g\" HaplotypeCaller -ERC GVCF --native-pair-hmm-threads 4 -I Clone$i.SORTED.MD.bam -O $VCF_DIR/Clone$i.g.vcf -R $REF \n",
" done"
]
},
{
"cell_type": "markdown",
"id": "3f006059-6271-45dc-a089-8c13b6bda946",
"metadata": {},
"source": [
"#### Check that the gvcf files have been created"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "aae87033-8412-4079-b1de-758dac08eb47",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"clone10.gvcf Clone13.g.vcf.idx Clone18.g.vcf density.csv\n",
"Clone10.g.vcf Clone14.g.vcf Clone18.g.vcf.idx genome.txt\n",
"clone10.gvcf.idx Clone14.g.vcf.idx Clone19.g.vcf output.vcf\n",
"Clone10.g.vcf.idx Clone15.g.vcf Clone19.g.vcf.idx output.vcf.idx\n",
"Clone11.g.vcf Clone15.g.vcf.idx Clone1.g.vcf rawSNP.vcf\n",
"Clone11.g.vcf.idx Clone16.g.vcf Clone1.g.vcf.idx rawSNP.vcf.idx\n",
"Clone12.g.vcf Clone16.g.vcf.idx Clone20.g.vcf\n",
"Clone12.g.vcf.idx Clone17.g.vcf Clone20.g.vcf.idx\n",
"Clone13.g.vcf Clone17.g.vcf.idx Clone2.g.vcf\n"
]
}
],
"source": [
"ls $VCF_DIR/"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0251d02d-1c88-4b9d-9f06-52f5cd2d841b",
"metadata": {
"tags": []
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "07a52085-ad7b-4833-b8c4-845085a2ff23",
"metadata": {
"tags": []
},
"source": [
"## __III - COMMON STEPS ANALYSES__ \n",
" \n",
"### Combine all the GVCF in on gvcf `gatk CombineGVCF` \n",
" \n",
"[See doc about CombineGVCF](https://gatk.broadinstitute.org/hc/en-us/articles/360037053272-CombineGVCFs)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "6afa0781-11db-4c1a-ad29-5d8c4dbe1efd",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" ################## J'ai généré la variable de type string : --variant Clone1.g.vcf --variant Clone2.g.vcf ############## \n",
"\n",
"\n",
"Using GATK jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar\n",
"Running:\n",
" java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar CombineGVCFs -R /home/jovyan/work/SV_DATA/REF/reference.fasta --variant Clone1.g.vcf --variant Clone2.g.vcf -O rawSNP.vcf\n",
"04:22:46.103 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_compression.so\n",
"04:22:46.307 INFO CombineGVCFs - ------------------------------------------------------------\n",
"04:22:46.308 INFO CombineGVCFs - The Genome Analysis Toolkit (GATK) v4.2.6.1\n",
"04:22:46.308 INFO CombineGVCFs - For support and documentation go to https://software.broadinstitute.org/gatk/\n",
"04:22:46.316 INFO CombineGVCFs - Initializing engine\n",
"04:22:46.780 INFO FeatureManager - Using codec VCFCodec to read file file:///home/jovyan/work/VCF/Clone1.g.vcf\n",
"04:22:46.807 INFO FeatureManager - Using codec VCFCodec to read file file:///home/jovyan/work/VCF/Clone2.g.vcf\n",
"04:22:46.862 INFO CombineGVCFs - Done initializing engine\n",
"04:22:46.889 INFO ProgressMeter - Starting traversal\n",
"04:22:46.890 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute\n",
"04:22:46.956 WARN ReferenceConfidenceVariantContextMerger - Detected invalid annotations: When trying to merge variant contexts at location Reference:39 the annotation MLEAC=[1, 0] was not a numerical value and was ignored\n",
"04:22:51.527 INFO ProgressMeter - Reference:995144 0.1 104470 1352070.8\n",
"04:22:51.527 INFO ProgressMeter - Traversal complete. Processed 104470 total variants in 0.1 minutes.\n",
"04:22:51.587 INFO CombineGVCFs - Shutting down engine\n",
"[June 22, 2022 4:22:51 AM UTC] org.broadinstitute.hellbender.tools.walkers.CombineGVCFs done. Elapsed time: 0.09 minutes.\n",
"Runtime.totalMemory()=1411907584\n"
]
}
],
"source": [
"BAM_ILL=\"/home/jovyan/work/BAM_ILL\"\n",
"VCF_DIR=\"/home/jovyan/work/VCF\"\n",
"REF=\"/home/jovyan/work/SV_DATA/REF/reference.fasta\"\n",
"\n",
"# change of work directory\n",
"cd $VCF_DIR\n",
"# Loop to inflate the --variant option\n",
"OPTION=\"\"\n",
"for i in {1..2}\n",
"do\n",
" OPTION=\"${OPTION} --variant Clone${i}.g.vcf\"\n",
"done\n",
"echo -e \"\\n ################## J'ai généré la variable de type string : $OPTION ############## \\n\\n\"\n",
"\n",
"# GATK\n",
"gatk CombineGVCFs -R $REF $OPTION -O rawSNP.vcf"
]
},
{
"cell_type": "markdown",
"id": "71d2f82d-dfc9-45ef-bca4-c186c356b4fe",
"metadata": {},
"source": [
"### Have a Look to the vcf created by combineGVCF "
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "9dbd61ee-0e7d-4bda-bc26-bfeb0830193f",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"##fileformat=VCFv4.2\n",
"##ALT=\n",
"##FILTER=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"Reference\t999881\t.\tA\t\t.\t.\tEND=999928\tGT:DP:GQ:MIN_DP:PL\t./.:9:27:9:0,27,330\t./.\n",
"Reference\t999929\t.\tC\t\t.\t.\tEND=999943\tGT:DP:GQ:MIN_DP:PL\t./.:9:24:8:0,24,326\t./.\n",
"Reference\t999944\t.\tT\t\t.\t.\tEND=999960\tGT:DP:GQ:MIN_DP:PL\t./.:8:21:7:0,21,285\t./.\n",
"Reference\t999961\t.\tC\t\t.\t.\tEND=999967\tGT:DP:GQ:MIN_DP:PL\t./.:7:18:7:0,18,270\t./.\n",
"Reference\t999968\t.\tT\t\t.\t.\tEND=999971\tGT:DP:GQ:MIN_DP:PL\t./.:7:15:6:0,15,225\t./.\n",
"Reference\t999972\t.\tT\t\t.\t.\t.\tGT:DP:GQ:MIN_DP:PL\t./.:6:12:6:0,12,180\t./.\n",
"Reference\t999973\t.\tG\t\t.\t.\tEND=999974\tGT:DP:GQ:MIN_DP:PL\t./.:6:9:6:0,9,135\t./.\n",
"Reference\t999975\t.\tT\t\t.\t.\tEND=999978\tGT:DP:GQ:MIN_DP:PL\t./.:6:6:5:0,6,90\t./.\n",
"Reference\t999979\t.\tT\t\t.\t.\tEND=999985\tGT:DP:GQ:MIN_DP:PL\t./.:4:3:2:0,3,45\t./.\n",
"Reference\t999986\t.\tC\t\t.\t.\tEND=1020001\tGT:DP:GQ:MIN_DP:PL\t./.:0:0:0:0,0,0\t./.\n"
]
}
],
"source": [
"head $VCF_DIR/rawSNP.vcf\n",
"tail $VCF_DIR/rawSNP.vcf"
]
},
{
"cell_type": "markdown",
"id": "e2ec071c-36de-498a-bc4b-2099bd1308e4",
"metadata": {},
"source": [
"### Compute the Genotypes - `gatk GenotypeGVCFs` \n",
"\n",
"GenotypeGVCFs uses the potential variants from the HaplotypeCaller and does the joint genotyping. It will look at the available information for each site from both variant and non-variant alleles across all samples, and will produce a VCF file containing only the sites that it found to be variant in at least one sample."
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "b6fbf4fe-d12f-4bf3-92b9-9aa446867e61",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using GATK jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar\n",
"Running:\n",
" java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xmx4g -jar /opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar GenotypeGVCFs -R /home/jovyan/work/SV_DATA/REF/reference.fasta -V /home/jovyan/work/VCF/rawSNP.vcf -O /home/jovyan/work/VCF/SNP.genotypeGVCF.vcf\n",
"04:25:33.845 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.2.6.1-1/gatk-package-4.2.6.1-local.jar!/com/intel/gkl/native/libgkl_compression.so\n",
"04:25:34.072 INFO GenotypeGVCFs - ------------------------------------------------------------\n",
"04:25:34.072 INFO GenotypeGVCFs - The Genome Analysis Toolkit (GATK) v4.2.6.1\n",
"04:25:34.072 INFO GenotypeGVCFs - For support and documentation go to https://software.broadinstitute.org/gatk/\n",
"04:25:34.081 INFO GenotypeGVCFs - Initializing engine\n",
"04:25:34.505 INFO FeatureManager - Using codec VCFCodec to read file file:///home/jovyan/work/VCF/rawSNP.vcf\n",
"04:25:34.544 INFO GenotypeGVCFs - Done initializing engine\n",
"04:25:34.596 INFO ProgressMeter - Starting traversal\n",
"04:25:34.596 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute\n",
"04:25:34.752 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position Reference:39 and possibly subsequent; at least 10 samples must have called genotypes\n",
"04:25:40.773 INFO ProgressMeter - Reference:993111 0.1 102607 996829.7\n",
"04:25:40.773 INFO ProgressMeter - Traversal complete. Processed 102607 total variants in 0.1 minutes.\n",
"04:25:40.780 INFO GenotypeGVCFs - Shutting down engine\n",
"[June 22, 2022 4:25:40 AM UTC] org.broadinstitute.hellbender.tools.walkers.GenotypeGVCFs done. Elapsed time: 0.12 minutes.\n",
"Runtime.totalMemory()=1387790336\n"
]
}
],
"source": [
"gatk --java-options \"-Xmx4g\" GenotypeGVCFs -R $REF -V $VCF_DIR/rawSNP.vcf -O $VCF_DIR/SNP.genotypeGVCF.vcf"
]
},
{
"cell_type": "markdown",
"id": "4629c532-f09c-4d1e-8e1d-9e58c78abc03",
"metadata": {},
"source": [
"### Have a Look to the vcf created by GenotypeGVCFs "
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "466b3284-5406-4701-b794-867e234c3c85",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"total 41244\n",
"-rw-r--r-- 1 jovyan users 8267535 Jun 21 22:36 Clone1.g.vcf\n",
"-rw-r--r-- 1 jovyan users 22384 Jun 21 22:36 Clone1.g.vcf.idx\n",
"-rw-r--r-- 1 jovyan users 2358535 Jun 21 22:39 Clone2.g.vcf\n",
"-rw-r--r-- 1 jovyan users 17736447 Jun 21 23:21 clone10.gvcf\n",
"-rw-r--r-- 1 jovyan users 37643 Jun 21 23:21 clone10.gvcf.idx\n",
"-rw-r--r-- 1 jovyan users 11068884 Jun 22 04:22 rawSNP.vcf\n",
"-rw-r--r-- 1 jovyan users 27643 Jun 22 04:22 rawSNP.vcf.idx\n",
"-rw-r--r-- 1 jovyan users 2684809 Jun 22 04:25 SNP.genotypeGVCF.vcf\n",
"-rw-r--r-- 1 jovyan users 4294 Jun 22 04:25 SNP.genotypeGVCF.vcf.idx\n",
"##fileformat=VCFv4.2\n",
"##ALT=\n",
"##FILTER=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"##FORMAT=\n",
"Reference\t998862\t.\tA\tC\t1901.06\t.\tAC=2;AF=1.00;AN=2;DP=46;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.67;SOR=0.874\tGT:AD:DP:GQ:PL\t1/1:0,46:46:99:1915,138,0\t./.:0,0\n",
"Reference\t998943\t.\tA\tT\t1760.06\t.\tAC=2;AF=1.00;AN=2;DP=46;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=33.97;SOR=0.784\tGT:AD:DP:GQ:PL\t1/1:0,44:44:99:1774,132,0\t./.:0,0\n",
"Reference\t999125\t.\tG\tT\t1711.06\t.\tAC=2;AF=1.00;AN=2;DP=45;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.01;SOR=0.892\tGT:AD:DP:GQ:PL\t1/1:0,42:42:99:1725,126,0\t./.:0,0\n",
"Reference\t999411\t.\tA\tC\t1423.06\t.\tAC=2;AF=1.00;AN=2;DP=36;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.58;SOR=0.869\tGT:AD:DP:GQ:PL\t1/1:0,35:35:99:1437,105,0\t./.:0,0\n",
"Reference\t999452\t.\tC\tA\t1253.06\t.\tAC=2;AF=1.00;AN=2;DP=32;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.80;SOR=1.214\tGT:AD:DP:GQ:PL\t1/1:0,31:31:93:1267,93,0\t./.:0,0\n",
"Reference\t999631\t.\tT\tC\t1141.06\t.\tAC=2;AF=1.00;AN=2;DP=30;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.71;SOR=0.906\tGT:AD:DP:GQ:PL\t1/1:0,29:29:87:1155,87,0\t./.:0,0\n",
"Reference\t999661\t.\tA\tC\t1078.06\t.\tAC=2;AF=1.00;AN=2;DP=27;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=27.67;SOR=1.302\tGT:AD:DP:GQ:PL\t1/1:0,27:27:81:1092,81,0\t./.:0,0\n",
"Reference\t999834\t.\tT\tA\t520.06\t.\tAC=2;AF=1.00;AN=2;DP=14;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=34.17;SOR=5.283\tGT:AD:DP:GQ:PL\t1/1:0,13:13:39:534,39,0\t./.:0,0\n",
"Reference\t999851\t.\tT\tC\t591.06\t.\tAC=2;AF=1.00;AN=2;DP=14;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.12;SOR=5.421\tGT:AD:DP:GQ:PL\t1/1:0,14:14:42:605,42,0\t./.:0,0\n",
"Reference\t999878\t.\tG\tT\t389.06\t.\tAC=2;AF=1.00;AN=2;DP=11;ExcessHet=0.0000;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=32.17;SOR=4.804\tGT:AD:DP:GQ:PL\t1/1:0,10:10:30:403,30,0\t./.:0,0\n"
]
}
],
"source": [
"ls -lrt\n",
"head $VCF_DIR/SNP.genotypeGVCF.vcf\n",
"tail $VCF_DIR/SNP.genotypeGVCF.vcf"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b4f88040-5973-4053-b672-53bda9d7f6d4",
"metadata": {},
"outputs": [],
"source": []
}
],
"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
}