added runs_per_week.sh and imporved compute_production_stats.py in order to plot per month yield. Updated docs to make Phil happy

vezzi · vezzi · commit 1952788fb4b1 · 2017-08-31T14:52:38.000+02:00
diff --git a/README.md b/README.md
@@ -23,36 +23,51 @@ Examples:
     - `python compute_undet_index_stats.py --config couch_db.yaml -- mode most_undet --instrument-type HiSeqX`
 
 
+### compute_undet_index_stats.py
+used to fetch stats about undermined indexes.
+This scripts queries statusdb x_flowcell_db  and fetch informaiton about runs.
+The following operations are supported:
 
-
-### DupRateTrends_from_charon.py
-Used to fetch stats from charon about duplication rate trends and number of sequenced human genomes
+ - check_undet_index: given a specific index checks all FCs and prints all FC and lanes where the indx appears as undetermined
+ - most_undet: outputs a summary about undetermiend indexes, printing the most 20 most occurring indexes for each instrument type
+ - single_sample_lanes: prints stats about HiSeqX lanes run with a single sample in it
+ - workset_undet: prints for each workset the FC, lanes and samples where the specified index has been found in undet. For each sample the plate position is printed.
+ - fetch_pooled_projects: returns pooled projects, that is projects that have been run in a pool. 
 
 #### Usage
-Example: `DupRateTrends_from_charon.py`
+Examples:
+ 
+  - compute for each workset the FC that contain a lane with index CTTGTAAT present in undet at least 0.5M times:
+    -  `python compute_undet_index_stats.py --config couch_db.yaml --index CTTGTAAT --mode workset_undet --min_occurences 500000` 
+ - Compute a list of the most occurring undetemriend indexes for HiSeqX runs:
+    - `python compute_undet_index_stats.py --config couch_db.yaml -- mode most_undet --instrument-type HiSeqX`
+
 
-```
-Usage: DupRateTrends_from_charon.py
 
-Options:
-  -h, --help            show this help message and exit
-  -t TOKEN, --token TOKEN
-                        Charon API Token. Will be read from the env variable
-                        CHARON_API_TOKEN if not provided
-  -u URL, --url URL     Charon base url. Will be read from the env variable
-                        CHARON_BASE_URL if not provided
-```
+
+
+### runs_per_week.sh
+Run on Irma prints a three columns:
+ 
+  - first column is the week number
+  - second column number of HiSeqX runs in that week
+  - seconf column number of HiSeq2500 runs in that week
+
+#### Usage
+Examp `runs_per_week.sh `
+
 
 
 ### compute_production_stats.py
-This scripts queries statusdb x_flowcelldb and project database and fetches informations about what organism have been sequenced. More in detail:
+This scripts queries statusdb x_flowcelldb and project database and fetches informations useful to plot trands and aggregated data. It can be run in three modalities:
+
+         - production-stats: for each instrument type it prints number of FCs, number of lanes, etc. It then prints a summary of all stats
+         - instrument-usage: for each instrument type and year it prints different run set-ups and samples run with that set-up
+         - year-stats: cumulative data production by month
 
-- reports total number of lanes sequenced per year
-- reports total number of Human lanes and of Non-Human lanes sequenced (divided per instrument)
-- other stats...
 
 ##### Usage
-Example: `compute_production_stats.py --config couchdb.yaml`
+Example: `compute_production_stats.py --config couchdb.yaml --mode year-stats`
 ```
 Usage: compute_production_stats.py --config couchdb.yam
 
diff --git a/compute_production_stats.py b/compute_production_stats.py
@@ -235,10 +235,8 @@ def instrument_usage():
                                         }
     flowcell_db = couch["x_flowcells"]
     project_sequenced = {}
-    FC = 0
+    instrument_runs_per_week = {}
     for fc_doc in flowcell_db:
-        if FC > 50:
-            continue
         if 'RunInfo' not in flowcell_db[fc_doc]:
             continue
         instrument = flowcell_db[fc_doc]["RunInfo"]['Instrument']
@@ -285,6 +283,13 @@ def instrument_usage():
         month =  int(flowcell_db[fc_doc]['RunInfo']['Date'][2:4])
         day   = int(flowcell_db[fc_doc]['RunInfo']['Date'][4:6])
         date_seq = datetime(year , month , day )
+        if instrument not in instrument_runs_per_week:
+            instrument_runs_per_week[instrument] = {}
+        date_entry = "{}_{}".format(year, month) #date_seq.isocalendar()[1]) # year plus week number
+        if date_entry not in instrument_runs_per_week[instrument]:
+            instrument_runs_per_week[instrument][date_entry] = 1
+        else:
+            instrument_runs_per_week[instrument][date_entry] += 1
         for lane in projects_in_lanes:
             for project in projects_in_lanes[lane]:
                 projects[project]['lanes'] += 1
@@ -293,11 +298,8 @@ def instrument_usage():
                         projects[project]['date'] = date_seq
                 else:
                     projects[project]['date'] = date_seq
-        #FC += 1
     flowcell_db = couch["flowcells"]
     for fc_doc in flowcell_db:
-        if FC > 100:
-            continue
         if 'RunInfo' not in flowcell_db[fc_doc]:
             continue
         if 'Date' not in flowcell_db[fc_doc]['RunInfo']:
@@ -307,17 +309,15 @@ def instrument_usage():
             print "run {} too old".format(flowcell_db[fc_doc]['RunInfo']['Id'])
             continue
         if 'Instrument' not in flowcell_db[fc_doc]["RunInfo"]:
-            import pdb
-            pdb.set_trace()
+            print "ERROR: Instrument not found in RunInfo: how is this possible?"
+            exit
 
         instrument = flowcell_db[fc_doc]["RunInfo"]['Instrument']
         if 'illumina' not in flowcell_db[fc_doc]:
             print "Not illumina field found in doc {}".format(fc_doc)
             continue
         if 'Demultiplex_Stats' not in  flowcell_db[fc_doc]['illumina']:
             print "Not Demultiplex_Stats field found in doc {}".format(fc_doc)
-            import pdb
-            pdb.set_trace()
             continue
         if 'Barcode_lane_statistics' not in flowcell_db[fc_doc]['illumina']['Demultiplex_Stats']:
             print "Not Barcode_lane_statistics field found in doc {}".format(fc_doc)
@@ -359,6 +359,14 @@ def instrument_usage():
         month =  int(flowcell_db[fc_doc]['RunInfo']['Date'][2:4])
         day   = int(flowcell_db[fc_doc]['RunInfo']['Date'][4:6])
         date_seq = datetime(year , month , day )
+        if instrument not in instrument_runs_per_week:
+            instrument_runs_per_week[instrument] = {}
+        date_entry = "{}_{}".format(year, month) #date_seq.isocalendar()[1]) # year plus week number
+        if date_entry not in instrument_runs_per_week[instrument]:
+            instrument_runs_per_week[instrument][date_entry] = 1
+        else:
+            instrument_runs_per_week[instrument][date_entry] += 1
+
         for lane in projects_in_lanes:
             for project in projects_in_lanes[lane]:
                 projects[project]['lanes'] += 1
@@ -367,9 +375,23 @@ def instrument_usage():
                         projects[project]['date'] = date_seq
                 else:
                     projects[project]['date'] = date_seq
-        #FC += 1
-
 
+    years = [2013, 2014, 2015, 2016, 2017]
+    sys.stdout.write('date,')
+    for instrument in sorted(instrument_runs_per_week):
+        sys.stdout.write('{},'.format(instrument))
+    sys.stdout.write('\n')
+    for year in years:
+        for week in xrange(1,13):
+            date_to_search = "{}_{}".format(year, week)
+            sys.stdout.write('{},'.format(date_to_search))
+            for instrument in sorted(instrument_runs_per_week):
+                if date_to_search in instrument_runs_per_week[instrument]:
+                    sys.stdout.write('{},'.format(instrument_runs_per_week[instrument][date_to_search]))
+                else:
+                    sys.stdout.write('0,')
+            sys.stdout.write('\n')
+    sys.stdout.write('\n')
     
     years = [2013, 2014, 2015, 2016, 2017]
     sequencers_year_setup = {}
@@ -409,6 +431,58 @@ def instrument_usage():
             sys.stdout.write('\n')
 
 
+def year_bp_production():
+    couch       = setupServer(CONFIG)
+    db_names = ['flowcells', 'x_flowcells']
+    flowcells   = {}
+    production_stats = {}
+    for db_name in db_names:
+        flowcell_db = couch[db_name]
+        for fc_doc in flowcell_db:
+            if 'RunInfo' not in flowcell_db[fc_doc]:
+                continue
+            if 'Flowcell' not in flowcell_db[fc_doc]['RunInfo']:
+                continue 
+            fc_name = flowcell_db[fc_doc]['RunInfo']['Flowcell']
+            if fc_name in flowcells:
+                continue
+            else:
+                flowcells[fc_name] = 0
+            year  = int(flowcell_db[fc_doc]['RunInfo']['Date'][0:2])
+            month = int(flowcell_db[fc_doc]['RunInfo']['Date'][2:4])
+            if year < 12:
+                continue
+            yield_MBases = 0
+            if 'illumina' not in flowcell_db[fc_doc]:
+                continue
+            if 'Demultiplex_Stats' not in flowcell_db[fc_doc]['illumina']:
+                continue
+            if db_name == "x_flowcells":
+                if 'Flowcell_stats' not in flowcell_db[fc_doc]['illumina']['Demultiplex_Stats']:
+                    continue
+                if 'Yield (MBases)' not in flowcell_db[fc_doc]['illumina']['Demultiplex_Stats']['Flowcell_stats']:
+                    continue
+                yield_MBases = int(flowcell_db[fc_doc]['illumina']['Demultiplex_Stats']['Flowcell_stats']['Yield (MBases)'].replace(',', ''))
+            else:
+                if 'Barcode_lane_statistics' not in  flowcell_db[fc_doc]['illumina']['Demultiplex_Stats']:
+                    continue
+                for sample  in flowcell_db[fc_doc]['illumina']['Demultiplex_Stats']['Barcode_lane_statistics']:
+                    yield_MBases +=  int(sample['Yield (Mbases)'].replace(',', ''))
+
+            if year not in production_stats:
+                production_stats[year] = month_production = [0]*12
+            production_stats[year][month-1] += yield_MBases
+    sys.stdout.write(',')
+    for year in sorted(production_stats):
+        sys.stdout.write('{},'.format(year))
+    sys.stdout.write('\n')
+    for month in range(0,12,1):
+        sys.stdout.write('{},'.format(month+1))
+        for year in sorted(production_stats):
+            sys.stdout.write('{},'.format(production_stats[year][month]))
+        sys.stdout.write('\n')
+    sys.stdout.write('\n')
+
 
     
 def main(args):
@@ -422,7 +496,10 @@ def main(args):
     
     if args.mode == 'instrument-usage':
         instrument_usage()
-    
+
+    if args.mode == 'year-stats':
+        year_bp_production()
+
     
 
 
@@ -434,9 +511,10 @@ def main(args):
     parser = argparse.ArgumentParser("""This scripts queries statusdb x_flowcelldb and project database and fetches informations about what organisms have been sequenced. It can be run in the following modes:
          - production-stats: for each instrument type it prints number of FCs, number of lanes, etc. It then prints a summary of all stats
          - instrument-usage: for each instrument type and year it prints different run set-ups and samples run with that set-up
+         - year-stats: cumulative data production by month
         """)
     parser.add_argument('--config', help="configuration file", type=str,  required=True)
-    parser.add_argument('--mode', help="define what action needs to be executed", type=str, required=True, choices=('production-stats', 'instrument-usage'))
+    parser.add_argument('--mode', help="define what action needs to be executed", type=str, required=True, choices=('production-stats', 'instrument-usage', 'year-stats'))
 
     args = parser.parse_args()
     main(args)
diff --git a/runs_per_week.sh b/runs_per_week.sh
@@ -0,0 +1,54 @@
+function runs_per_week()
+{
+    local fc_root_folders=/proj/ngi2016003/*/
+    local week=$1 year=$2
+    local week_num_of_Jan_1 week_day_of_Jan_1
+    local first_Mon
+    local date_fmt="+%y %m %d"
+    local date_ftm_year="+%y"
+    local date_fmt_month="+%m"
+    local date_fmt_day="+%d"
+
+    week_num_of_Jan_1=$(date -d $year-01-01 +%W)
+    week_day_of_Jan_1=$(date -d $year-01-01 +%u)
+
+    if ((week_num_of_Jan_1)); then
+        first_Mon=$year-01-01
+    else
+        first_Mon=$year-01-$((01 + (7 - week_day_of_Jan_1 + 1) ))
+    fi
+
+
+    YEAR=${year: -2}
+    MONTH_START=$(date -d "$first_Mon +$((week - 1)) week" "$date_fmt_month")
+    MONTH_END=$(date -d "$first_Mon +$((week - 1)) week + 6 day" "$date_fmt_month")
+    DAY_START=$(date -d "$first_Mon +$((week - 1)) week" "$date_fmt_day")
+    DAY_END=$(date -d "$first_Mon +$((week - 1)) week + 6 day" "$date_fmt_day")
+
+    mon=$(date -d "$first_Mon +$((week - 1)) week" "$date_fmt")
+    sun=$(date -d "$first_Mon +$((week - 1)) week + 6 day" "$date_fmt")
+    DAYS=()
+    if [ $MONTH_START -ne $MONTH_END ];
+    then
+        DAYS=($(seq -f "$YEAR$MONTH_START%02g" $DAY_START 1 31) $(seq -f "$YEAR$MONTH_END%02g" 1 1 $DAY_END))
+    else
+        DAYS=($(seq -f "$YEAR$MONTH_START%02g" $DAY_START 1 $DAY_END))
+    fi
+
+    RUNS_PER_WEEK_X=0
+    RUNS_PER_WEEK_nonX=0
+    for DAY in "${DAYS[@]}" ; do
+        DAY_RUNS_X=`ls -d $fc_root_folders/$DAY*_ST* 2> null | wc -l`
+        RUNS_PER_WEEK_X=`expr $DAY_RUNS_X + $RUNS_PER_WEEK_X`
+        DAY_RUNS_nonX=`ls -d $fc_root_folders/$DAY* 2> null | grep -v ST | grep -v 000000  | wc -l`
+        RUNS_PER_WEEK_nonX=`expr $DAY_RUNS_nonX + $RUNS_PER_WEEK_nonX`
+    done
+    echo $week $RUNS_PER_WEEK_X $RUNS_PER_WEEK_nonX
+
+}
+
+CURRENT_YEAR=`date +"%Y"`
+CURRENT_WEEK=`date +"%V"`
+for WEEK in `seq 1 1 $CURRENT_WEEK`; do
+        runs_per_week $WEEK $CURRENT_YEAR
+done
