CLOUD-BASED SOMATIC VARIANT CALLING ON CD56+ CTC DATA USING AWS EC2 & S3

1. Introduction

This internship project was undertaken to gain hands-on experience with cloud-based bioinformatics infrastructure and reproducible pipeline development. The primary goals were to learn how to provision and manage AWS EC2 instances and set up S3 buckets for scalable data storage and result management.

As the applied use case, a somatic variant calling pipeline was executed on publicly available NGS data from CD56+ Circulating Tumour Cells (CTCs) of a Small Cell Lung Cancer (SCLC) patient, sourced from Ricordel et al. (2023). The sample ERR6473446 was used as the input dataset.

The GATK Best Practices pipeline, covering alignment, pre-processing, Mutect2 variant calling, and SnpEff annotation, was replicated on an AWS EC2 instance following the methodology described in the reference study. The biological findings of the original paper were used as a qualitative benchmark to verify that the pipeline was correctly implemented, rather than to generate new scientific conclusions.

2. Objective

The objective of this internship project was to learn and implement a reproducible, cloud-based somatic variant calling workflow on public NGS data using:

• AWS EC2 for compute.
• AWS S3 for storage.
• A standard GATK-based pipeline for variant calling, following the configuration used in Ricordel et al. (2023).

The primary focus was on infrastructure setup and pipeline execution, not on generating novel biological interpretations.

3. Cloud Infrastructure Setup

3.1 AWS Account and IAM

Created an AWS account and configured IAM user with least-privilege access (permissions for EC2, S3). Generated an access key and secret key for CLI use and configured them with aws configure.

3.2 S3 Bucket

Created an S3 bucket to store: 1) Raw FASTQ files (optional, if mirrored from ENA). 2) Intermediate BAM/VCF outputs. 3) Final annotated VCFs and result summaries.

Example Bucket name: s3://adilsukumar/

Folder structure:

• data/ - (To store raw FASTQ files, filtered & unfiltered VCF files, etc.)
• results/ - (Intermediate BAM/VCF outputs and final annotated summaries)
• reference/ - (Human reference genomes (GRCh38))

3.3 EC2 Instance Configuration

To handle the computational demands of NGS alignment and variant calling, a virtual server was deployed and configured as follows:

1. Instance Launch and Access

• OS: Ubuntu Linux.
• Instance Deployment: The virtual server was launched via the AWS Management Console.
• Permissions Management: To ensure the private key remained secure and accessible only by the owner, permissions were adjusted from the local terminal using:
  

  chmod 400 [path_to_key].pem

• Remote Login: A secure connection was established from the local terminal to the cloud server using SSH (Secure Shell):
  

  ssh -i [path_to_key_pemfile] ubuntu@[public_ip]

2. Environment and Software Setup

• System Update & Base Dependencies: A standard sudo apt update was performed to ensure all package repositories and system libraries were current. Essential system utilities including bash, curl, and wget were verified to facilitate data downloading and script execution.
• Environment Management: Mamba: To ensure reproducible dependency management and avoid library conflicts, Mamba (a fast C++ implementation of Conda) was installed. A dedicated bioinformatics environment was created to host the primary pipeline tools.
• Manual Bioinformatics Tool Installation: Core preprocessing, alignment, and postprocessing tools were installed directly on the Ubuntu environment or via Mamba. This allowed for direct execution and specific versioning of the following software:
  • SRA Toolkit (prefetch, fasterq-dump): Used to retrieve raw datasets from NCBI and convert SRA files into FASTQ format.
  • FastQC & MultiQC: Utilized to check the quality of individual samples and aggregate those metrics into a single run report.
  • Fastp: Employed for ultra-fast all-in-one FASTQ preprocessing and quality control (QC), generating reports in HTML format.
  • BWA (Burrows-Wheeler Aligner): Used for aligning the sequence data to the GRCh38 human reference genome.
  • Samtools: Utilized for converting SAM files to BAM format, as well as sorting and indexing alignment files.
  • GATK (Genome Analysis Toolkit): The primary engine for variant calling, used for: MarkDuplicates (identifying PCR duplicates), BaseRecalibrator / ApplyBQSR (systemic error correction) and Mutect2 (somatic variant calling).
  • SnpEff: A genetic variant annotation and functional effect prediction tool. It was used to categorize the effects of variants (such as deletions or frame-shift mutations) directly within the Linux environment.
• Verification: Software functionality was confirmed for all manually installed tools using the --version flag (e.g., samtools --version) to ensure they were correctly added to the system path and ready for the pipeline.

3.4 Storage Synchronization

Efficient data handling is critical in NGS workflows due to the massive size of raw FASTQ, intermediate BAM, and final VCF files.

• S3 Bucket Linking: The EC2 working directory was linked to an Amazon S3 bucket, creating a seamless bridge between the local high-speed processing environment and persistent cloud storage.
• Connectivity and Backups: Connectivity was verified using aws s3 ls. To prevent data loss upon instance termination, the aws s3 sync command was utilized to transfer processed reports and genomic variants directly to the cloud repository.

4. Methodology

A standard bioinformatics pipeline was utilized to process the raw sequencing data within a Linux environment. The workflow proceeded as follows:

4.1 Commands Used in Linux

Step 1: Data Download

Downloaded raw reads (SRA/FASTQ) using prefetch and fasterq-dump of the sample and of the Human Reference Genome.

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR647/ERR6473446/ERR6473446_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR647/ERR6473446/ERR6473446_2.fastq.gz

wget https://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
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

Step 2: Quality Control

FastQC was run to assess the quality of the raw reads. The analysis indicated high sequence quality (Phred scores > 30); therefore, no trimming was required.

fastqc ERR6473446_1.fastq.gz ERR6473446_2.fastq.gz

Step 3: Alignment to the Reference Genome

The high-quality reads were aligned to the human reference genome (GRCh38) using BWA-MEM, a standard tool for the accurate alignment of short-read sequencing data.

bwa mem -t 16 ../reference/Homo_sapiens_assembly38.fasta \
ERR6473446_1.fastq.gz ERR6473446_2.fastq.gz > aligned.sam

Step 4: Conversion, Sorting & Indexing

The SAM file was converted to the binary BAM format using samtools. The resulting BAM file was sorted by genomic coordinates and indexed.

samtools view -bS aligned.sam > aligned.bam
samtools sort aligned.bam -o aligned.sorted.bam
samtools index aligned.sorted.bam

Step 5: Adding Read Group Information

Read group metadata, essential for GATK-based analysis, was added to the sorted BAM file.

gatk AddOrReplaceReadGroups \
  -I aligned.sorted.bam \
  -O aligned.rg.bam \
  -RGID 1 \
  -RGLB lib1 \
  -RGPL ILLUMINA \
  -RGPU unit1 \
  -RGSM TUMOR_SAMPLE

Step 6: Marking Duplicate Reads

PCR and sequencing duplicates were flagged using GATK MarkDuplicates to prevent bias in variant calling.

gatk MarkDuplicates \
  -I aligned.rg.bam \
  -O aligned.dedup.bam \
  -M duplicates.metrics.txt
samtools index aligned.dedup.bam

Step 7: Base Quality Score Recalibration (BQSR)

BQSR was performed in two sub-steps to correct systematic errors introduced during the sequencing process.

7.1 Creating the Recalibration Table

gatk BaseRecalibrator \
  -I aligned.dedup.bam \
  -R ~/project/reference/Homo_sapiens_assembly38.fasta \
  --known-sites ~/project/reference/Homo_sapiens_assembly38.dbsnp138.vcf \
  -O recal_data.table

7.2 Applying the Recalibration

gatk ApplyBQSR \
  -R /home/ubuntu/project/reference/Homo_sapiens_assembly38.fasta \
  -I /home/ubuntu/project/results/aligned.dedup.bam \
  --bqsr-recal-file /home/ubuntu/project/results/recal_data.table \
  -O /home/ubuntu/project/results/recal.bam

Step 8: Somatic Variant Calling

The recalibrated BAM file (recal.bam) served as the definitive input. Somatic variants were identified using GATK Mutect2, with gnomAD serving as the germline resource.

gatk --java-options "-Xmx16g" Mutect2 \
  -R ~/project/reference/Homo_sapiens_assembly38.fasta \
  -I ~/project/results/recal.bam \
  -tumor TUMOR_SAMPLE \
  --germline-resource ~/project/reference/af-onlygnomad.hg38.vcf.gz \
  --native-pair-hmm-threads 20 \
  -O ~/project/results/tumor.unfiltered.vcf.gz

Step 9: Filtering Variant Calls

Mutect2 output was filtered to remove low-confidence and artifact calls before annotation.

gatk FilterMutectCalls \
  -R ~/project/reference/Homo_sapiens_assembly38.fasta \
  -V ~/project/results/tumor.unfiltered.vcf.gz \
  -O ~/project/results/tumor.filtered.vcf.gz

Step 10: Annotation

SnpEff was used to annotate the Mutect2 generated VCF file. This analysis pinpointed the affected genes and forecasted how the variants might alter proteins, for instance, by truncating them or swapping one amino acid for another.

java -Xmx8g -jar snpEff.jar \
  -v GRCh38.99 \
  ~/project/results/tumor.unfiltered.vcf.gz > ~/project/results/tumor.ann.vcf

Step 11: Filtering for High & Moderate Impact Variants

SnpSift was used to filter the annotated VCF for variants with significant predicted biological impact (High or Moderate), removing low-impact or modifier variants.

java -Xmx4g -jar SnpSift.jar filter \
  "(ANN[0].IMPACT = 'HIGH') | (ANN[0].IMPACT = 'MODERATE')" \
  ~/project/results/tumor.ann.vcf > ~/project/results/tumor.high_mod.vcf

Step 12: Extracting Coding Variants

Further isolation of variants specifically altering the protein-coding sequence, such as missense mutations, stop gained (nonsense), and frameshift variants.

java -Xmx4g -jar SnpSift.jar filter \
  "ANN[0].EFFECT has 'missense_variant' \
  | ANN[0].EFFECT has 'stop_gained' \
  | ANN[0].EFFECT has 'stop_lost' \
  | ANN[0].EFFECT has 'start_lost' \
  | ANN[0].EFFECT has 'frameshift_variant'" \
  ~/project/results/tumor.ann.vcf > ~/project/results/tumor.coding_effect.vcf

Step 13: Syncing Results to S3

This makes the run reproducible and sharable without keeping the EC2 instance alive.

aws s3 sync /project/results s3://adilsukumar/results/

4.2 Tools & Commands Used

This analysis was conducted on an AWS EC2 instance (Ubuntu) for containerized tool management and S3 for scalable data storage.

5. Results

The genomic analysis of the CD56+ CTC sample (ERR6473446) revealed a highly mutated and heterogeneous tumor profile.

5.1 High Tumor Mutational Burden (TMB): CTCs vs. Biopsy

A critical finding from the reference study and our analysis is the massive discrepancy in tumor mutation burden between liquid biopsy (CTCs) and solid tissue biopsy. CTCs capture the "sum" of all metastases, whereas biopsies only capture a small fraction.

5.2 SnpEff Analysis: Variant Consequences

Using SnpEff, we classified the identified somatic mutations based on their predicted functional impact. The results demonstrate a high concentration of high-impact structural alterations in the CD56+ CTC sample. Specifically, 85.7% of key driver mutations are Indel/Frameshift alterations and 14.3% are Nonsense (Stop Gained) mutations. These disruptive mutations are hallmark characteristics of the aggressive genomic instability found in SCLC.

• Indel/Frameshift Mutations (85.7%): The dominant mutation type, affecting TP53, CREBBP, EP300, COL22A1, PTEN, and NOTCH3. These lead to truncated proteins or loss-of-function, driving uncontrolled survival and metastatic potential.
• Non-Sense Mutations (14.3%): Specifically identified in RB1, introducing a premature stop codon that abolishes cell-cycle brakes, triggering aggressive tumor proliferation.
• Missense Mutations (0% in Key Driver Subset): While common in the broader genome, the primary drivers are defined by more disruptive structural alterations.

5.3 Mutational Signatures: The "Smoking Gun" vs. Evolution

The analysis revealed distinct mutational processes driving the tumor at different sites. The Tissue Biopsy a massive spike in C>A mutations (the direct "smoking gun" of tobacco damage), while the CTCs show a shift toward C>T mutations, indicating that as the cells circulate, they acquire new evolutionary mutations unrelated to the original smoking damage.

5.4 Key Driver Gene Alterations

6. Discussion

The implementation of this NGS pipeline on CD56+ CTCs (Sample ERR6473446) provides a critical "liquid window" into the evolving SCLC genome, capturing insights that traditional solid biopsies often miss due to necrosis or limited sampling.

6.1 Capturing Tumor Heterogeneity

The pipeline results demonstrate that a single tissue biopsy fails to capture the tumor’s full genetic diversity. The analyzed CD56+ CTCs revealed a Tumor Mutational Burden (TMB) significantly higher (20 to 30 times) than matched tissue biopsies. By circulating through the vasculature, these CTCs effectively represent the "sum" of multiple metastatic sites rather than a localized, single-site fraction.

6.2 Tracking Genomic Evolution

The mutational signature analysis indicates a clear shift in mutational processes. While the primary tissue biopsy is dominated by C>A mutations, the CTCs show a transition toward C>T mutations. This suggests that as the cancer progresses and cells enter circulation, they acquire new evolutionary mutations independent of the original environmental triggers.

6.3 Genomic "Priming" for Metastasis

The pipeline identified a specific subset of high-impact driver mutations in genes such as COL22A1 (extracellular matrix destabilization), NOTCH3 (disruption of neuroendocrine differentiation), and PTEN (loss of PI3K pathway regulation). The presence of these alterations suggests that these circulating cells are genomically "primed" for survival in the bloodstream and subsequent metastatic seeding.

7. Conclusion

This project successfully established a custom bioinformatics pipeline to process CD56+ Liquid Biopsy data for Small Cell Lung Cancer. By leveraging AWS infrastructure and GATK Best Practices, the workflow bypassed the limitations of necrotic tissue samples to provide a comprehensive genomic profile.

The pipeline validated the obligate loss of SCLC’s defining tumor suppressors, accurately detecting:

• TP53: Indel and frameshift mutations that nullify DNA repair.
• RB1: Nonsense mutations that abolish cell-cycle control.
• Epigenetic & Signaling Regulators: High-impact mutations in CREBBP, EP300, and the PTEN pathway.

Ultimately, this computational framework demonstrates that CD56+ CTC analysis is a superior method for tracking SCLC evolution, offering a scalable approach for future personalized medicine and real-time tumor monitoring.

References

1. Ricordel et al. (2023). Genomic characteristics and clinical significance of CD56+ circulating tumor cells in small cell lung cancer. Scientific Reports. PMC9984363. Sample Accession ERR6473446.
2. Broad Institute. GATK Best Practices: Somatic Short Mutation Calling.
3. European Nucleotide Archive (ENA). Sample Accession ERR6473446.