<jupyter_start><jupyter_text>A deep dive into TOB-WGS' pipelineChapter 1: what does `ReblockGVCF` actually do ? Let's examine sample TOB1520 (which is an outlier for some QC metrics)<jupyter_code>import hail as hl;

# All datasets in TOB-WGS are using GRCh38
hl.init(default_reference='GRCh38');<jupyter_output>Running on Apache Spark version 3.1.2
SparkUI available at http://loic-notebook.australia-southeast1-a.c.notebooks-314505.internal:4040
Welcome to
     __  __     <>__
    / /_/ /__  __/ /
   / __  / _ `/ / /
  /_/ /_/\_,_/_/_/   version 0.2.67-bafea6b18247
LOGGING: writing to /home/jupyter/tob-deepdive/hail-20210618-0309-0.2.67-bafea6b18247.log<jupyter_text>Import a GVCF and explore structure and content<jupyter_code>gvcf = hl.import_vcf('gs://cpg-tob-wgs-test/gvcf/batch1/TOB1520.g.vcf.gz', min_partitions=12, force_bgz=True)
gvcf.describe()
gvcf.info.END.show()
hl.summarize_variants(gvcf)<jupyter_output>2021-06-18 03:11:32 Hail: INFO: Coerced sorted dataset<jupyter_text>Worth noting:of 35,550,357 variants, 30,075,685 are homozygous reference blocksLet's load the GVCF after `ReblockGVCF` has been performed<jupyter_code>rb_gvcf = hl.import_vcf('gs://cpg-tob-wgs-test-tmp/joint-calling/v2/hail/batch/728b87/1/output_gvcf.g.vcf.gz', force_bgz=True, min_partitions=12)
rb_gvcf.describe();
hl.summarize_variants(rb_gvcf);<jupyter_output>----------------------------------------
Global fields:
    None
----------------------------------------
Column fields:
    's': str
----------------------------------------
Row fields:
    'locus': locus<GRCh38>
    'alleles': array<str>
    'rsid': str
    'qual': float64
    'filters': set<str>
    'info': struct {
        AC: array<int32>, 
        AF: array<float64>, 
        AN: int32, 
        AS_BaseQRankSum: array<float64>, 
        AS_FS: array<float64>, 
        AS_InbreedingCoeff: array<float64>, 
        AS_MQ: array<float64>, 
        AS_MQRankSum: array<float64>, 
        AS_QD: array<float64>, 
        AS_QUALapprox: str, 
        AS_RAW_BaseQRankSum: str, 
        AS_RAW_MQ: str, 
        AS_RAW_MQRankSum: str, 
        AS_RAW_ReadPosRankSum: str, 
        AS_ReadPosRankSum: array<float64>, 
        AS_SB_TABLE: str, 
        AS_SOR: array<float64>, 
        AS_VarDP: str, 
        BaseQRankSum: float64, 
        DP: int32, 
        DS: bool, 
        END: int32, 
   [...]<jupyter_text>These are the annotations added by ReblockGVCF        AC, AF, AN: ???? why would we need that on a single sample ?        AS_BaseQRankSum: array,         AS_FS: array,        AS_MQ: array,         AS_MQRankSum: array,        AS_QD: array,         AS_QUALapprox: str,        AS_ReadPosRankSum: array,        AS_SOR: array,         AS_VarDP: str,        FS: float64,         MQ: float64,          MQ_DP: int32,         QD: float64,         QUALapprox: int32,         RAW_GT_COUNT: array,         SOR: float64,         VarDP: int32        and the annotations removed by ReblockGVCF        MLEAC: array,         MLEAF: array Additionally we lost more than 5e6 variants: not only homref blocks but also SNPs, etc.Let's have a look at what we lost:<jupyter_code>lost = gvcf.anti_join_rows(rb_gvcf.rows())
hl.summarize_variants(lost)<jupyter_output>2021-06-18 03:23:08 Hail: INFO: Coerced sorted dataset
2021-06-18 03:23:57 Hail: INFO: Coerced sorted dataset<jupyter_text>We lost mostly homref blocks, as expected, but also other kind of variants, such as 425,204 SNPsLet's inspect some of these lost variants (more than 2 alleles means it is not a homref block)<jupyter_code>lost.filter_rows(hl.len(lost.alleles) == 3).show(50,1,truncate=30)<jupyter_output>2021-06-18 03:30:59 Hail: INFO: Coerced sorted dataset
2021-06-18 03:31:56 Hail: INFO: Coerced sorted dataset
2021-06-18 03:32:46 Hail: INFO: Coerced sorted dataset<jupyter_text>Why did we lose for example chr1:16571 ? a 0/0 with GQ=54. Did it get merged in an homref block ?Note: the phased GT does not match the GT ? how is that possible ? (example: chr1:109262)Note: all these variant are homref. Did we lose any non homref variants ?<jupyter_code>lost = lost.annotate_rows(nonhomref=hl.agg.count_where(lost.GT.is_non_ref())>0)
lost.aggregate_rows(hl.agg.count_where(lost.nonhomref))<jupyter_output>2021-06-18 03:42:11 Hail: INFO: Coerced sorted dataset
2021-06-18 03:43:02 Hail: INFO: Coerced sorted dataset<jupyter_text>156,380 non homref variants were lost, let's examine some of them:<jupyter_code>lost.filter_rows(lost.nonhomref).show(50,1, width=10)<jupyter_output>2021-06-18 04:05:42 Hail: INFO: Coerced sorted dataset
2021-06-18 04:06:41 Hail: INFO: Coerced sorted dataset
2021-06-18 04:07:33 Hail: INFO: Coerced sorted dataset
StarcoderdataJupyterStructuredCleanDedup