vcf_into_sqlite: make WITHOUT ROWID optional

Author: mlin

loaders/vcf_into_sqlite.cc

@@ -45,7 +45,8 @@ vector<map<string, string>> extract_hrecs(bcf_hdr_t *hdr, const char *key,
 }
 
 string schemaDDL(const string &table_prefix, vector<map<string, string>> &info_hrecs,
-                 vector<map<string, string>> &format_hrecs, int ploidy) {
+                 vector<map<string, string>> &format_hrecs, int ploidy,
+                 bool genotypes_without_rowid) {
 
     OStringStream out;
     out << "CREATE TABLE " << table_prefix
@@ -120,7 +121,10 @@ string schemaDDL(const string &table_prefix, vector<map<string, string>> &info_h
                 }
             }
         }
-        out << "\n, PRIMARY KEY (variant_rowid, sample_id)) WITHOUT ROWID";
+        out << "\n, PRIMARY KEY (variant_rowid, sample_id))";
+        if (genotypes_without_rowid) {
+            out << " WITHOUT ROWID";
+        }
     }
 
     return string(out.Get());
@@ -612,30 +616,33 @@ void help() {
         << '\n'
         << "vcf_into_sqlite [options] in.vcf|- out.db" << '\n'
         << "Options: " << '\n'
-        << "  --table-prefix PREFIX  prefix to the name of each table created" << '\n'
-        << "  --ploidy N             set max ploidy => # GT columns (default 2)" << '\n'
-        << "  --no-gri               skip genomic range indexing" << '\n'
-        << "  -l,--level LEVEL       database compression level (-7 to 22, default 6)" << '\n'
-        << "  -q,--quiet             suppress progress information on standard error" << '\n'
-        << "  -h,--help              show this help message" << '\n'
+        << "  --table-prefix PREFIX      prefix to the name of each table created" << '\n'
+        << "  --ploidy N                 set max ploidy => # GT columns (default 2)" << '\n'
+        << "  --genotypes-without-rowid  make the genotypes table WITHOUT ROWID (advantageous if there are few/no FORMAT fields)"
+        << '\n'
+        << "  --no-gri                   skip genomic range indexing" << '\n'
+        << "  -l,--level LEVEL           database compression level (-7 to 22, default 6)" << '\n'
+        << "  -q,--quiet                 suppress progress information on standard error" << '\n'
+        << "  -h,--help                  show this help message" << '\n'
         << '\n';
 }
 
 int main(int argc, char *argv[]) {
     string table_prefix, infilename, outfilename;
-    bool gri = true, progress = true;
+    bool gri = true, progress = true, genotypes_without_rowid = false;
     int level = 6, ploidy = 2;
 
     static struct option long_options[] = {{"help", no_argument, 0, 'h'},
                                            {"quiet", no_argument, 0, 'q'},
                                            {"level", required_argument, 0, 'l'},
                                            {"ploidy", required_argument, 0, 'p'},
+                                           {"genotypes-without-rowid", no_argument, 0, 'R'},
                                            {"table-prefix", required_argument, 0, 't'},
                                            {"no-gri", no_argument, 0, 'G'},
                                            {0, 0, 0, 0}};
 
     int c;
-    while (-1 != (c = getopt_long(argc, argv, "hqGl:t:p:", long_options, nullptr))) {
+    while (-1 != (c = getopt_long(argc, argv, "hqGRl:t:p:", long_options, nullptr))) {
         switch (c) {
         case 'h':
             help();
@@ -646,6 +653,9 @@ int main(int argc, char *argv[]) {
         case 'G':
             gri = false;
             break;
+        case 'R':
+            genotypes_without_rowid = true;
+            break;
         case 't':
             table_prefix = optarg;
             break;
@@ -751,7 +761,8 @@ int main(int argc, char *argv[]) {
 
         // formulate & apply DDL
         // TODO: allow --append
-        string ddl = schemaDDL(table_prefix, info_hrecs, format_hrecs, ploidy);
+        string ddl =
+            schemaDDL(table_prefix, info_hrecs, format_hrecs, ploidy, genotypes_without_rowid);
         progress &&cerr << ddl << endl;
         db->exec(ddl);
 

test/genomicsqlite_big_tests.wdl

@@ -1,6 +1,6 @@
 version 1.0
 
-workflow genomicsqlite_integration_tests {
+workflow genomicsqlite_big_tests {
     input {
         String git_revision = "main"
 
@@ -74,6 +74,7 @@ task build {
 task test_sam {
     input {
         File reads
+        String cram_ref_fa_gz = "ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz"
 
         File genomicsqlite_py
         File libgenomicsqlite_so
@@ -84,8 +85,9 @@ task test_sam {
 
     command <<<
         set -euxo pipefail
+        TMPDIR=${TMPDIR:-/tmp}
         apt-get -qq update
-        DEBIAN_FRONTEND=noninteractive apt-get -qq install -y sqlite3 samtools tabix libzstd1
+        DEBIAN_FRONTEND=noninteractive apt-get -qq install -y sqlite3 samtools tabix libzstd1 pigz wget
 
         cp ~{genomicsqlite_py} /usr/lib/python3.8/genomicsqlite.py
         cp ~{libgenomicsqlite_so} /usr/local/lib/libgenomicsqlite.so
@@ -94,8 +96,19 @@ task test_sam {
         cp ~{sam_into_sqlite} /usr/local/bin/sam_into_sqlite
         chmod +x /usr/local/bin/sam_into_sqlite
 
+        reads_file='~{reads}'
+        if [[ $reads_file == *.cram ]]; then
+            # CRAM given; make BAM to take reference downloads out of the timings
+            wget -nv -O - '~{cram_ref_fa_gz}' | pigz -dc > "${TMPDIR}/cram_ref.fa"
+            samtools faidx "${TMPDIR}/cram_ref.fa"
+            bam_file="${TMPDIR}/$(basename reads_file .cram).bam"
+            time samtools view -T "${TMPDIR}/cram_ref.fa" -h -O BAM -@ 8 "$reads_file" > "$bam_file"
+            reads_file=$bam_file
+        fi
+        >&2 ls -l "$reads_file"
+
         # load database
-        time sam_into_sqlite "~{reads}" "~{dbname}"
+        time sam_into_sqlite "$reads_file" "~{dbname}"
         >&2 ls -l "~{dbname}"
 
         # GRI query
@@ -168,7 +181,7 @@ task test_vcf {
         chr = genomicsqlite.get_reference_sequences_by_name(dbconn)
         query = 'SELECT count(*) FROM ' + genomicsqlite.genomic_range_rowids_sql(dbconn, 'variants')
         print(query, file=sys.stderr)
-        row = next(dbconn.execute(query, (chr['chr12'].rid,111803912,111804012))) #(chr['chr21'].rid, 34787801, 35049344)))
+        row = next(dbconn.execute(query, (chr['chr21'].rid, 34787801, 35049344))) #(chr['chr12'].rid,111803912,111804012)))
         print(f'result = {row[0]}', file=sys.stderr)
         print(row[0])
         EOF