The Analytic Garden

 On November 22, 2022, I downloaded a metadata file containing 13,932,236 records from GISAID  in order to view the growth of several emerging SARS-CoV-2 variants of interest. First, the results. I plotted the growth of the variants BA.2.75, BQ.1, BQ.1.1, BQ.1.18, CL.1, and XBB in the USA from June 1 to the present.  BQ.1. and BQ.1.1 are growing in the data. It's unclear yet whether they will have a large impact. We'll see after the Thanksgiving holidays. Memory Issues The metadata file is large, 9.8 GB on the disk. Loading into R on Windows takes a little over two minutes and used about 8 GB for the dataframe. For the purposes of demonstration, I'm using Windows 10 22H2, RStudio 2022.07.2, and R 4.2.2. > system.time(meta_data <- read_tsv( 'data/metadata.tsv' , name_repair = 'universal' )) user system elapsed 197.84 9.44 146.04 > object.size(meta_data) 8087686552 bytes > RStudio uses 9.6 GB once the data is loaded. PS...

 The world doesn't really need another description of how to code logistic regression. A good description of how to implement logistic regression can already be found here . In addition, there are many great packages for logistic regression in Python,  sklearn.linear_model. LogisticRegression ; R glm ; Julia GLM ; and many more.  I started following this course on Udemy. The course began with a brief discussion of logistic regression. I have used logistic regression techniques many times, but I didn't have a clear idea of how to implement it. I thought I might as well try. What follows is a very simple implementation of binary logistic regression . Logistic Regression - Who lives, who dies? Consider the following data. It's from the R package alr4 . It describes the fate of the infamous Donner party . It consists of 91 observations of five variables. The important columns for our purposed are age , y (survival), and sex . We want to know how did age and sex affect s...

 Here are the updates for the five county NY Capital District area. I don't trust the case data. Cases are under reported, but here is what is being reported. Cases are trending up. This is worrisome because I have said I and others believe cases are underreported. Home tests are not being reported and mild cases probably aren't even being tested. Hospital data may be more reliable.  Hospitalization seem low, but the counties vary in population.  Albany and Columbia counties have relatively high rates. Luckily deaths are low. Code for downloading data and producing plots is available on  GitHub .

 President Biden says the pandemic is over . Maybe, but still ~400 Americans are dying each day from COVID-19 and excess deaths are still 10% above “expected.” It seems SARS-CoV-2 is not done with us. See the wonderful  Your Local Epidemiologist blog for details. In fact, I would like to put in a plug for Katelyn Jetelina's blog . It's one of the best sources for explanations of the current epidemiological situations that I have encountered. On September 13 2022, I downloaded a metadata file from  GISAID  containing 13,061,086 rows and 22 columns.  For this post, I take a brief look at the progress of the BA variants in the US. It's hard to know if the GISAID data gives a representative sample of the COVID-19 situation as, like other dats sources, it's likely that the actual situation is underreported. Still, the data you have is the dats that you have. Here's what the data shows for variants of concern. As we have seen before the Omicron variant is firm...

This post describes a little experiment in loading large file with Julia, Python, and R. On September 13, 2022, I downloaded a metadata file from GISAID . The file 13,061,086 rows and 22 columns. As a file for the PC, it's big. My usual approach to analyzing GISAID data is using a combination of R for counting and plotting variables such as linages over time and using Python, particularly BioPython, for manipulating sequence data.  When using R, I tend to use the tidyverse  tools for manipulating tabular data. Transforming dataframes by piping then through functions seems like a natural approach. Julia, the other hand, is often fast enough that writing simple loops to manipulate data is feasible and can lead to simpler more readable code. Loading a file with more than 13 million rows is slow in R. I wondered if Julia or Python/Pandas could do better.  What follows is an unscientific exercise in reading a large tab delimited file. All the tests were run on a PC with a...

 For a recent project I needed sequence regions upstream (preceding then 5' end of the gene) of a set of orthologous genes. The orthologs for a gene of interest are obtained from  https://www.ncbi.nlm.nih.gov/gene . For example, searching for  JAK2 orthologs  at that site yields a table of JAK2 genes for a large number of species. After selecting species, the ortholog table can be downloaded. Fetching Genomes Since I wanted to analyze a number of different genes, I decided to automate the process of getting the upstream regions. The first step was to fetch the GenBank records for the genomes of the selected species. The GenBank IDs for each species are included in the downloaded ortholog table. Fetching genomes is straightforward, if a bit slow. It uses Pandas to read the ortholog Table from NCBI and BioPython.Entrez to download the complete GenBank record for the genome. def main (): args = GetArgs() genome_path = args . genome_path ortholog_tab...