update the coverage calculation method

Author: tayabsoomro

frontend/src/modules/analysis/analysis-data/analysis-data.component.tsx

@@ -10,9 +10,12 @@ const AnalysisDataComponent: FunctionComponent<IAnalysisDataProps> = ({ data })
     const [coverageData, setCoverageData] = useState<any[]>([]);
     const [listenerRunning, setListenerRunning] = useState(false);
     const [error, setError] = useState("");
-    const [metric, setMetric] = useState<'avg_depth' | 'breadth'>('avg_depth');
+    const [metric, setMetric] = useState<'fold_coverage'>('fold_coverage');
     const [showConfirmModal, setShowConfirmModal] = useState(false);
 
+    // Extract threshold from the first query, default to 100 if invalid
+    const threshold = data.data.queries[0]?.threshold ? parseFloat(data.data.queries[0].threshold) : 100;
+
     useEffect(() => {
         socket.emit('check_fastq_file_listener', { projectId: analysisData.data.projectId });
 
@@ -72,29 +75,27 @@ const AnalysisDataComponent: FunctionComponent<IAnalysisDataProps> = ({ data })
     const formatCoverageData = () => {
         const refs = [...new Set(coverageData.map(d => d.reference))];
         const times = [...new Set(coverageData.map(d => d.timestamp))].sort();
-        const startTime = new Date(times[0]).getTime(); // Earliest timestamp in milliseconds
-    
-        // Define the header with column types and roles
-        const header = [{ type: 'number', label: 'Elapsed Time (s)' }];
+        const startTime = new Date(times[0]).getTime();
+
+        const header = [{ type: 'number', label: 'Elapsed Time (s)', role: '' }];
         refs.forEach(ref => {
-            header.push({ type: 'number', label: ref });
-            header.push({ type: 'string', label: '' }); // Tooltip column, no 'role' property
+            header.push({ type: 'number', label: ref, role: '' });
+            header.push({ type: 'string', label: 'for', role: 'tooltip' });
         });
-    
-        // Create data rows with elapsed time and tooltips
+
         const rows = times.map(time => {
-            const elapsedSeconds = (new Date(time).getTime() - startTime) / 1000; // Convert to seconds
-            const row: any[] = [elapsedSeconds];
+            const elapsedSeconds = (new Date(time).getTime() - startTime) / 1000;
+            const row: (number | string)[] = [elapsedSeconds];
             refs.forEach(ref => {
                 const entry = coverageData.find(d => d.timestamp === time && d.reference === ref);
-                const y = entry ? entry[metric] : 0;
-                const tooltip = `Time: ${time}\n${ref}: ${y.toFixed(2)}`;
+                const y = entry ? entry.fold_coverage : 0;
+                const tooltip = `Time: ${time}\n${ref}: ${y.toFixed(2)}x`;
                 row.push(y);
                 row.push(tooltip);
             });
             return row;
         });
-    
+
         return [header, ...rows];
     };
 
@@ -136,6 +137,23 @@ const AnalysisDataComponent: FunctionComponent<IAnalysisDataProps> = ({ data })
         setShowConfirmModal(false);
     };
 
+    // Compute number of references for series indexing
+    const refs = [...new Set(coverageData.map(d => d.reference))];
+    const numRefs = refs.length;
+
+    const chartOptions = {
+        title: 'Fold Coverage Over Time',
+        hAxis: { title: 'Elapsed Time (s)' },
+        vAxis: { title: 'Fold Coverage (x)', minValue: 0 },
+        legend: { position: 'bottom' },
+        colors: ['#00B0BD', '#004E5A', '#FF6A45', '#27AE60'],
+        chartArea: { width: '80%', height: '70%' },
+        animation: { startup: true, duration: 1000, easing: 'out' },
+        series: !isNaN(threshold) ? {
+            [numRefs]: { lineDashStyle: [4, 4], color: 'red', lineWidth: 2 }
+        } : {}
+    };
+
     return (
         <div className="nano-analysis-container">
             <h2 className="nano-section-title">Analysis Dashboard</h2>
@@ -188,10 +206,10 @@ const AnalysisDataComponent: FunctionComponent<IAnalysisDataProps> = ({ data })
                     <h3 className="nano-card-title">Coverage Over Time</h3>
                     <select
                         value={metric}
-                        onChange={(e) => setMetric(e.target.value as 'avg_depth' | 'breadth')}
+                        onChange={(e) => setMetric(e.target.value as 'fold_coverage')}
                         className="form-control w-auto d-inline-block ml-2"
                     >
-                        <option value="avg_depth">Average Depth</option>
+                        <option value="fold_coverage">Average Depth</option>
                         <option value="breadth">Breadth of Coverage (%)</option>
                     </select>
                 </div>
@@ -200,19 +218,7 @@ const AnalysisDataComponent: FunctionComponent<IAnalysisDataProps> = ({ data })
                         <Chart
                             chartType="LineChart"
                             data={formatCoverageData()}
-                            options={{
-                                title: metric === 'avg_depth' ? 'Average Coverage Depth Over Time' : 'Breadth of Coverage Over Time',
-                                hAxis: { title: 'Elapsed Time (s)' },
-                                vAxis: { 
-                                    title: metric === 'avg_depth' ? 'Average Depth (reads/position)' : 'Breadth (%)', 
-                                    minValue: 0,
-                                    maxValue: metric === 'breadth' ? 100 : undefined
-                                },
-                                legend: { position: 'bottom' },
-                                colors: ['#00B0BD', '#004E5A', '#FF6A45', '#27AE60'],
-                                chartArea: { width: '80%', height: '70%' },
-                                animation: { startup: true, duration: 1000, easing: 'out' }
-                            }}
+                            options={chartOptions}
                             width="100%"
                             height="400px"
                         />

frontend/src/modules/setup/setup-steps/database-setup/additional-sequences-setup/additional-sequences-setup.component.tsx

@@ -5,7 +5,7 @@ type IKeys = "name" | "file" | "threshold" | "alert";
 
 const AdditionalSequencesSetupComponent: FunctionComponent<IDatabaseSetupConstituent> = ({ updateConfig }) => {
     const [queries, setQueries] = useState([
-        { name: "", file: "", threshold: "", current_breadth: 0, alert: false }
+        { name: "", file: "", threshold: "", current_fold_change: 0, alert: false, header: "" }
     ]);
 
     const handleDataChange = (idx: number, key: IKeys) => (evt: React.ChangeEvent<HTMLInputElement>) => {
@@ -23,7 +23,7 @@ const AdditionalSequencesSetupComponent: FunctionComponent<IDatabaseSetupConstit
             alert("Please fill all fields in the current query before adding a new one.");
             return;
         }
-        setQueries((prev) => [...prev, { name: "", file: "", threshold: "", current_breadth: 0, alert: false }]);
+        setQueries((prev) => [...prev, { name: "", file: "", threshold: "", current_fold_change: 0, alert: false, header: "" }]);
     };
 
     const handleRemoveQuery = (idx: number) => () => {
@@ -87,11 +87,10 @@ const AdditionalSequencesSetupComponent: FunctionComponent<IDatabaseSetupConstit
                                     id="thresholdText"
                                     name="thresholdText"
                                     type="number"
-                                    placeholder="Breadth Coverage % Threshold"
+                                    placeholder="Fold Coverage Threshold (x)"
                                     onChange={handleDataChange(i, "threshold")}
                                     className="form-control"
                                     min="0"
-                                    max="100"
                                 />
                             </div>
                             <div className="col-sm-1">

frontend/src/modules/setup/setup-steps/database-setup/additional-sequences-setup/additional-sequences-setup.interfaces.ts

@@ -6,8 +6,9 @@ type IQuery = {
     name: string,
     file: string,
     threshold: string,
-    current_breadth: number,
-    alert: boolean
+    current_fold_change: number,
+    alert: boolean,
+    header: string,
 }
 
 export type {

frontend/src/modules/setup/setup.component.tsx

@@ -16,7 +16,7 @@ import {IAlertNotifSetupInput} from "./setup-steps/alert-notif-setup/alert-notif
 
 const qrs: IAdditionalSequences = {
     queries: [
-        {name: "", file: "", threshold: "", current_breadth: 0, alert: false}
+        {name: "", file: "", threshold: "", current_fold_change: 0, alert: false, header: ""},
     ]
 };
 

server/app/main/routes.py

@@ -255,13 +255,12 @@ def get_coverage():
             lines = f.readlines()[1:]  # Skip header
         data = []
         for line in lines:
-            timestamp, ref, avg_depth, breadth, read_count = line.strip().split(',')
+            timestamp, ref, fold_coverage, read_count = line.strip().split(',')
             name = ref_to_name.get(ref, ref)  # Map reference to alert sequence name
             data.append({
                 'timestamp': timestamp,
                 'reference': name,
-                'avg_depth': float(avg_depth),
-                'breadth': float(breadth),
+                'fold_coverage': float(fold_coverage),
                 'read_count': int(read_count)
             })
         return jsonify(data)

server/app/main/utils/FileHandler.py

@@ -32,7 +32,7 @@ def __init__(self, app_loc: str):
         self.file_type = self.config.get('fileType', 'FASTQ')
         if not os.path.exists(self.coverage_file):
             with open(self.coverage_file, 'w') as f:
-                f.write("timestamp,reference,avg_depth,breadth,read_count\n")
+                f.write("timestamp,reference,fold_coverage,read_count\n")
 
     def on_moved(self, event):
         self.on_any_event(event)
@@ -150,7 +150,7 @@ def merge_bam(self, new_bam: str):
             logger.error(f"Error sorting/indexing merged BAM: {e}")
 
     def calculate_and_record_coverage(self):
-        """Calculate and record average depth, breadth of coverage, and read count per reference."""
+        """Calculate and record fold coverage (average depth) and read count per reference."""
         try:
             bam = pysam.AlignmentFile(self.merged_bam, "rb", check_sq=False)
             if not bam.has_index():
@@ -159,22 +159,23 @@ def calculate_and_record_coverage(self):
             coverage_data = {}
             for ref in bam.references:
                 ref_length = bam.lengths[bam.references.index(ref)]
-                pileup = bam.pileup(ref)
-                depths = [p.nsegments for p in pileup if p.nsegments > 0]
-                covered_positions = len(depths)
-                avg_depth = sum(depths) / len(depths) if depths else 0
-                breadth = (covered_positions / ref_length) * 100 if ref_length > 0 else 0
+                # Get coverage depth across all positions for this reference
+                coverage = bam.count_coverage(ref)
+                # Sum depths across all bases (A, C, G, T) at each position
+                total_depth = sum(sum(cov) for cov in coverage)
+                # Calculate fold coverage as total depth divided by reference length
+                fold_coverage = total_depth / ref_length if ref_length > 0 else 0
                 read_count = bam.count(ref)  # Count reads mapping to this reference
-                coverage_data[ref] = {"avg_depth": avg_depth, "breadth": breadth, "read_count": read_count}
-                # Check if breadth exceeds the threshold and trigger alert if necessary
-                print(f"Reference: {ref}, Avg Depth: {avg_depth:.2f}, Breadth: {breadth:.2f}%, Read Count: {read_count}")
-                self.check_breadth_alert(ref, breadth)
+                coverage_data[ref] = {"fold_coverage": fold_coverage, "read_count": read_count}
+                print(f"Reference: {ref}, Fold Coverage: {fold_coverage:.2f}x, Read Count: {read_count}")
+                # Update alert check to use fold coverage if needed
+                self.check_fold_coverage_alert(ref, fold_coverage)
             bam.close()
 
             timestamp = time.strftime("%Y-%m-%d %H:%M:%S", time.gmtime())
             with open(self.coverage_file, 'a') as f:
                 for ref, cov in coverage_data.items():
-                    f.write(f"{timestamp},{ref},{cov['avg_depth']},{cov['breadth']},{cov['read_count']}\n")
+                    f.write(f"{timestamp},{ref},{cov['fold_coverage']},{cov['read_count']}\n")
             logger.debug(f"Coverage and read counts recorded at {timestamp}")
 
             # Emit coverage update via Socket.IO
@@ -186,49 +187,22 @@ def calculate_and_record_coverage(self):
         except Exception as e:
             logger.error(f"Error calculating coverage: {e}")
 
-    def check_breadth_alert(self, ref: str, breadth: float):
-        """Check if breadth exceeds threshold and send alerts."""
-        print(f"Checking breadth alert for {ref} with breadth {breadth:.2f}%")
-        alertinfo_cfg_file = os.path.join(self.app_loc, 'alertinfo.cfg')
-        try:
-            with open(alertinfo_cfg_file, 'r') as f:
-                alertinfo_cfg_data = json.load(f)
-            queries = alertinfo_cfg_data.get("queries", [])
-            device = alertinfo_cfg_data.get("device", "")
-            print(alertinfo_cfg_data)
-            for query in queries:
-                print(f"Checking query: {query}")
-                if ref == query.get("header", ""):
-                    print(f"Matched query: {query}")
-                    threshold = float(query.get("threshold", 0))
-                    current_breadth = query.get("current_breadth", 0)
-                    # Check alert condition before updating current_breadth
-                    if breadth >= threshold and current_breadth < threshold:
-                        alert_str = f"Alert: {query['name']} breadth coverage reached {breadth:.2f}% (threshold: {threshold}%)"
-                        logger.critical(alert_str)
-                        if device:
-                            LinuxNotification.send_notification(device, alert_str)
-                        # Send alerts if configured
-                        if "alertNotifConfig" in alertinfo_cfg_data:
-                            print("Sending email")
-                            email_config = alertinfo_cfg_data['alertNotifConfig']
-                            subject = "nanoCAS Alert"
-                            body = alert_str
-                            send_email(subject, body, email_config)
-                        # Send SMS via Twilio if configured in .env
-                        if os.getenv('TWILIO_ACCOUNT_SID'):
-                            print("Sending SMS")
-                            send_sms(alert_str)
-                        else:
-                            print("Twilio not configured; skipping SMS")
-                            logger.debug("Twilio not configured in .env; SMS skipped")
-                    # Update current_breadth after checking
-                    if breadth > current_breadth:
-                        query["current_breadth"] = breadth
-                        with open(alertinfo_cfg_file, 'w') as f:
-                            json.dump(alertinfo_cfg_data, f, indent=4)
-                else:
-                    print(f"No match for query: {query}")
-                    logger.debug(f"No match for query: {query}")
-        except Exception as e:
-            logger.error(f"Error checking breadth alert: {e}")
+    def check_fold_coverage_alert(self, ref: str, fold_coverage: float):
+        """Check if fold coverage exceeds threshold and send alerts."""
+        with open(os.path.join(self.app_loc, 'alertinfo.cfg'), 'r') as f:
+            alertinfo_cfg_data = json.load(f)
+        queries = alertinfo_cfg_data.get("queries", [])
+        device = alertinfo_cfg_data.get("device", "")
+        for query in queries:
+            if ref == query.get("header", ""):
+                threshold = float(query.get("threshold", 0))
+                if fold_coverage >= threshold:
+                    alert_str = f"Alert: {query['name']} fold coverage reached {fold_coverage:.2f}x (threshold: {threshold}x)"
+                    logger.critical(alert_str)
+                    if device:
+                        LinuxNotification.send_notification(device, alert_str)
+                    if "alertNotifConfig" in alertinfo_cfg_data:
+                        email_config = alertinfo_cfg_data['alertNotifConfig']
+                        send_email("nanoCAS Alert", alert_str, email_config)
+                    if os.getenv('TWILIO_ACCOUNT_SID'):
+                        send_sms(alert_str)