iPSC QC Pipeline

This pipeline includes:

• Processing raw fastq files (bulk RNA sequencing) using fastqc, STAR alignment, and RSEM
• GATK germline variant calling
• Cancer mutation calling (in select oncogenes)
• eSNPKaryotyping
• Mycoplasma detection (Mycoplasma fermentans and Mycoplasma orale and Mycoplasma hyorhinis)
• PACNet classification
• Outlier assessment

This pipeline expects a directory of fastq.gz (paired end: R1 and R2) files from a bulk RNA sequencing run.

See wiki for more details.

---

Table of Contents

• Running the Pipeline
• Requirements
• Inputs
• Outputs
• Directory Structure
• Notes
• Sources and Citations

---

Running the Pipeline

Example Command

#Download docker image as singularity image (for HPC compatibility):

singularity pull ipsc_qc.sif docker://kjakubiak16/ipsc_qc:latest

#Run the singularity command. First bind the directories, then supply bound directories as arguments.

singularity exec \
  -B ./path/to/fastqs:/data \
  -B ./path/to/ref_dir:/ref \
  -B ./path/to/COSMIC_DB:/cosmic \
  -B ./path/to/output_dir:/output \
  ipsc_qc.sif \
  /pipeline/pipeline_runner.sh \
    --fastq_dir /data \
    --output_dir /output \
    --ref_dir /ref \
    --cosmic_dir /cosmic
    --keep_files true

Requirements

Reference COSMIC Database Files (provide in --cosmic_dir argument; default is ./cosmic)

Must be provided by user

• Cosmic_CancerGeneCensus_Tsv_v101_GRCh38.tar
  Cosmic Database; MUST BE DOWNLOADED BY USER AND PROVIDED IN COSMIC_DIR ARGUMENT. If not, cancer mutation calling will not be performed.

Reference Files (place in ref/ directory; default is ./ref)

• chr/ Directory with dbSNP build 142 common SNP files (GTF format), one per chromosome. Human: X = 23, Y = 24. Source: UCSC Table Browser, snp142common table.
• genes.txt File containing a list of oncogenes for assessment. Lives in this repo; download into your /ref dir (the one you provide to the --ref_dir argument), and add whichever genes you want analyzed against COSMIC.
• GTEX and RSEM References (Homo_sapiens_assembly38_noALT_noHLA_noDecoy.fasta, Homo_sapiens_assembly38_noALT_noHLA_noDecoy.fasta.fai, Homo_sapiens_assembly38_noALT_noHLA_noDecoy.dict, gencode.v39.GRCh38.genes.collapsed_only
  (requester pays enabled; cannot be universally downloaded within pipeline))

If not included in user-provided reference directory (ref_dir), the following files will be automatically downloaded to ref_dir:

• star_index_oh75 reference files for STAR
• rsem_reference refDict
• Homo_sapiens_assembly38.dbsnp138.vcf dbSnpVcf
• Homo_sapiens_assembly38.dbsnp138.vcf.idx dbSnpVcfIndex
• Mills_and_1000G_gold_standard.indels.hg38.vcf.gz knownVcf 1
• Homo_sapiens_assembly38.known_indels.vcf.gz knownVcf 2
• Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi knownVcfIndex 1
• Homo_sapiens_assembly38.known_indels.vcf.gz.tbi knownVcfIndex 2
• gencode.v39.GRCh38.genes.collapsed_only.gtf annotationsGTF
• GCF_000420105.1_ASM42010v1_genomic.fna.gz
  Mycoplasma orale genome; downloaded from wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/420/105/GCF_000420105.1_ASM42010v1/GCF_000420105.1_ASM42010v1_genomic.fna.gz
• GCF_003704055.1_ASM370405v1_genomic.fna.gz Mycoplasma fermentans genome; downloaded from wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/704/055/GCF_003704055.1_ASM370405v1/GCF_003704055.1_ASM370405v1_genomic.fna.gz
• Hs_expTrain_Jun-20-2017.rda
  PACNet training data; downloaded from aws s3 cp s3://cellnet-rnaseq/ref/cnproc/HS/Hs_expTrain_Jun-20-2017.rda . --no-sign-request
• Hs_stTrain_Jun-20-2017.rda
  PACNet training data; downloaded from aws s3 cp s3://cellnet-rnaseq/ref/cnproc/HS/Hs_stTrain_Jun-20-2017.rda . --no-sign-request
• GCF_000001405.26/GCF_000001405.26_GRCh38_genomic.fna
  NCBI reference genome FASTA file; downloaded from ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.26_GRCh38.p13

Tools

The pipeline depends on the following command line tools:

Must be installed by user (Ensure these tools are installed or available via a module system (ex. HPC users)):

• singularity
• parallel
• squashfuse

Inputs

Argument	Description	Default
--ref_dir	Directory for reference files	/refs
--fastq_dir	Directory containing input files (fastq.gz sequencing files) for a run	./fastq
--output_dir	Path to desired output directory	./outputs
--cosmic_dir	Directory containing downloaded COSMIC database files.	cosmic
--project	Project name	default_project
--keep_files	Keep bulky intermediate files (VCFs, BAMs, etc.)?	true

Note:

• File base names should be unique for each sample and consistent across file types (e.g., sample1.vcf, sample1.bam, sample1.genes.results, sample1_R1.fastq.gz, sample1_R2.fastq.gz)
• Pipeline will accept fastq files with base names formatted as SAMPLE1_01 & SAMPLE1_02 or SAMPLE1_R1 & SAMPLE1_R2 or SAMPLE1_<string>_<string>_R1_<string> & SAMPLE1_<string>_<string>_R2_<string> with file formats .fastq.gz or .fq.gz. See table below for more details. If your fastq files do not follow these formats, please reformat your files.

Acceptable fastq file formats:

Example File Name	Notes
samplename_1.fastq.gz	Simple numeric read identifier
samplename_2.fq.gz	Supports both .fastq.gz and .fq.gz
samplename_R1.fastq.gz	Standard Illumina R1/R2 notation
samplename_R2_001.fastq.gz	Trailing _001 ignored
sample-R1_001.fq.gz	Hyphen as separator supported
samplename-S7_L002_R2_001.fastq.gz	Handles _S and _L lane info
24246R-15-15_S0_L001_R1_001.fastq.gz	Complex, mixed characters, still parsed correctly

Key parsing rules:

• Sample name is everything before optional _S or _L and before R1/R2 or 1/2.
• Read number is extracted from R1, R2, 1, or 2 (case-insensitive).
• Supports .fastq.gz and .fq.gz.
• Extra numbering like _001 is ignored.
• Hyphens or underscores as separators are handled.

---

Outputs

Output Directory Structure

output_dir
  ├── logs
  │    ├── ekaryo_sample_01.log
  │    ├── ekaryo_sample_02.log
  │    ├── cancer_mutation_mapping_sample_01.log
  │    ├── cancer_mutation_mapping_sample_02.log
  │    ├── html_summary_sample_01.log
  │    ├── html_summary_sample_02.log
  │    ├── pacnet_all.log
  │    ├── outliers_sample_01.log
  │    ├── outliers_sample_02.log
  │    ├── sample_01_wdl2.log
  │    ├── sample_02_wdl2.log
  │    ├── sample_01_bulk_preprocess.log
  │    └── sample_02_bulk_preprocess.log  
  ├── pacnet
  │    ├── PACNet_heatmap.png
  │    ├── classifier.rda
  │    ├── query_meta.csv
  │    ├── query_matrix.csv
  │    ├── metrics.json
  │    ├── classification_validation_hm.png
  │    └── classification_scores.csv
  ├── outlier_analysis
  │    └── PCA_pacnet_scores.pdf
  ├── sample_01
  │    ├── mycoplasma
  │    │   ├── sample_01_flagstat.txt
  │    │   ├── sample_01_sorted.bam.bai
  │    │   ├── sample_01_sorted.bam
  │    │   ├── sample_01_bowtie2.log
  │    │   ├── mycoplasma_alignment_summary.png
  │    │   ├── mycoplasma_alignment_summary.pdf
  │    │   └── mycoplasma_alignment_stats.tsv
  │    ├── eSNPKaryotyping
  │    │   ├── sample_01_PlotGenome.png
  │    │   └── sample_01_variantTable.csv
  │    ├── cosmic_calling
  │    │   ├── CancerMutationPlot.pdf
  │    │   └── sample_01_CancerMutations.tsv
  │    ├── variant_calling
  │    │   ├── star_out
  │    │   │   ├── sample_01.Aligned.toTranscriptome.out.bam
  │    │   │   ├── sample_01.Aligned.out.bam
  │    │   │   ├── sample_01.Aligned.sortedByCoord.out.bam
  │    │   │   ├── sample_01.Aligned.sortedByCoord.out.bam.bai
  │    │   │   ├── sample_01.ReadsPerGene.out.tab.gz
  │    │   │   ├── sample_01.Chimeric.out.sorted.bam
  │    │   │   └── sample_01.Chimeric.out.sorted.bam.bai
  │    │   ├── RSEM_outputs
  │    │   │   ├── sample_01.rsem.isoforms.results.gz
  │    │   │   ├── sample_01.rsem.genes.results.gz
  │    │   │   └── sample_01.rsem.stat  
  │    │   ├── QC_outputs
  │    │   │   ├── sample_01_R1.fastqc_data.txt.gz
  │    │   │   ├── sample_01_R2.fastqc_data.txt.gz
  │    │   │   ├── sample_01_R1_fastqc.zip
  │    │   │   ├── sample_01_R2_fastqc.zip
  │    │   │   ├── sample_01_R1_fastqc.html
  │    │   │   └── sample_01_R2_fastqc.html
  │    │   ├── fastqc_out
  │    │   │   ├── sample_01_R1.fastqc_data.txt.gz
  │    │   │   ├── sample_01_R2.fastqc_data.txt.gz
  │    │   │   ├── sample_01_R1_fastqc.zip
  │    │   │   ├── sample_01_R2_fastqc.zip
  │    │   │   ├── sample_01_R1_fastqc.html
  │    │   │   └── sample_01_R2_fastqc.html
  │    │   └── Mark_duplicates_outputs
  │    │   │   ├── sample_01.Aligned.sortedByCoord.out.md.bam.bai
  │    │   │   ├── sample_01.marked_dup_metrics.txt
  │    │   │   └── sample_01.Aligned.sortedByCoord.out.md.bam
  │    ├── preprocessing
  │    ├── logs
  │    ├── variant_calling
  │    │   ├── sample_01.variant_filtered.vcf.gz
  │    │   ├── sample_01.variant_filtered.vcf.gz.tbi
  │    │   ├── sample_01.aligned.duplicates_marked.recalibrated.bam
  │    │   └── sample_01.aligned.duplicates_marked.recalibrated.bam.bai
  └── sample_02

Directory structure

Basic pipeline repo structure

src
    └── basic_pipeline
        ├── docker_image_structure
        ├── Dockerfile_ipsc_qc
        ├── modules
        │   ├── cancer_mutation_calling
        │   │   └── COSMIC_cancer_mutation_calling.r
        │   ├── eSNPKaryotyping
        │   │   └── run_eSNPKaryotyping.R
        │   ├── mycoplasma_detection
        │   │   └── detect_mycoplasma.sh
        │   ├── outlier_detection
        │   │   └── outlier_detection.R
        │   ├── PACNet
        │   │   └── run_pacnet.R
        │   ├── preprocessing
        │   │   └── wdlplay
        │   │       ├── LICENSE.md
        │   │       ├── Makefile
        │   │       ├── MANIFEST.in
        │   │       ├── README.md
        │   │       ├── requirements.txt
        │   │       ├── settings.json
        │   │       ├── setup.cfg
        │   │       ├── setup.py
        │   │       ├── versioneer.py
        │   │       ├── warp-pipelines
        │   │       │   ├── bulk_RNAseq_preprocess
        │   │       │   │   ├── RNAseq_pipeline_fastq.wdl
        │   │       │   │   ├── run_wdl.py
        │   │       │   │   └── wdl
        │   │       │   │       ├── fastqc.wdl
        │   │       │   │       ├── markduplicates.wdl
        │   │       │   │       ├── rnaseqc2.wdl
        │   │       │   │       ├── rsem.wdl
        │   │       │   │       └── star.wdl
        │   │       │   └── GATK_variant_calling
        │   │       │       ├── gatk4-rna-best-practices.wdl
        │   │       │       ├── gatk4-rna-germline-calling_run.py
        │   │       │       └── gatk4-rna-germline-variant-calling.inputs.json
        │   │       └── wdlplay
        │   │           ├── __init__.py
        │   │           ├── _version.py
        │   │           ├── db.json
        │   │           └── wdlplayer.py
        │   └── summary_writer.R
        ├── pipeline_runner.sh
        ├── ref_files
        │   └── genes.txt
        ├── report_builder.R
        └── tarballs
            └── CellNet_master.tar.gz

---

Notes

• This pipeline is memory-intensive. Implementation on HPC is recommended.
• To save space, set the --keep_files tag to false, especially if you do not intend to use the aligned BAM files or VCF files for other analyses.
• WDL scripts will utilize cromwell, but a caper server will not be started. Backend is local. All packages (including caper) and scripts are contained within the image; user does not need to install anything beyond what is included in the requirements section.

Interpretation Guidelines:

Based solely on bulk RNA sequencing data, this pipeline should show you comprehensive QC metrics for each iPSC line that you run. The interpretation of several of these metrics relies heavily on comparison with the other lines in your batch, and should give the user the ability to decide which of their lines to use for downstream experiments and analyses. The main purpose of the pipeline is not to provide complex or definitive predictions of differentiation outcomes, but to provide the user with enough metrics to make an informed decision on which lines to move forward with. Additionally, it provides a baseline framework for assessing the quality of iPSC lines, with the intention of standardizing which metrics are used for selection across labs. This will make the study of iPSC models and their efficacy more diagnosable and failures less subjective across studies.

• The PACNet module will let the user know if their cell line has deviated from the pluripotent state, or has begun to express lineage-specific markers. If one or multiple lines show an abnormal score in the “esc” category, this should be apparent in the PACNet heatmap as well as the summary heatmap. If a line has a lower score than the other lines in this category, it may no longer be pluripotent, therefore not suited for differentiation.

• The cancer mutation module detects mutations in cancer-associated genes–these are not uncommon and may oftentimes be benign, but if one or multiple lines have substantially higher numbers of mutations in genes like TP53, EGFR, and BRCA1, then this is easily identifiable from the summary heatmap. If a line has an abundance of mutations in these genes, it carries a higher likelihood of becoming cancerous.

• The eSNPKaryotyping module will detect chromosomal abnormalities, and the fraction abnormal karyotype metric shown in the summary heatmap will inform the user how much of the line’s genome has been affected by these abnormalities. If the fraction of abnormal SNPs in a given line is higher than the others, as seen in the summary heatmap, this may signal genomic instability, unpredictable differentiation outcomes, and/or tumorigenic potential.

• The mycoplasma detection module will assess the percent alignment of a line’s genome to the genomes of three species of mycoplasma. If a line has a substantial amount of alignment with any of the mycoplasma genomes–as depicted in the summary heatmap–then it is likely contaminated and therefore not suited for further analyses.

• The outlier analyses–based on PCA of counts and of esc scores–provide a way of determining which lines have significantly deviated from the others in terms of overall gene expression and expression of pluripotent markers. Lines shown to be outliers by this analytical framework may not produce optimal results, compared to other lines in the batch.

• The summary heatmap shows the metrics described above all at once; if a cell line is compromised, it may show aberrant scores across multiple modules, thereby alerting the user that it may not be suitable for downstream experiments and analyses. However, the user should use their discretion in assessing which modules produce aberrant scores for their lines, and which metrics may not be useful for their specific purposes.

---

Citations

• COSMIC: the Catalogue Of Somatic Mutations In Cancer: Tate JG, Bamford S, Jubb HC, Sondka Z, Beare DM, Bindal N, Boutselakis H, Cole CG, Creatore C, Dawson E, Fish P, Harsha B, Hathaway C, Jupe SC, Kok CY, Noble K, Ponting L, Ramshaw CC, Rye CE, Speedy HE, Stefancsik R, Thompson SL, Wang S, Ward S, Campbell PJ, Forbes SA. COSMIC: the Catalogue Of Somatic Mutations In Cancer. Nucleic Acids Res. 2019 Jan 8;47(D1):D941-D947. doi: 10.1093/nar/gky1015. PMID: 30371878; PMCID: PMC6323903.
• PACNet (Platform-Agnostic CellNet): Lo EKW, Velazquez JJ, Peng D, Kwon C, Ebrahimkhani MR, Cahan P. Platform-agnostic CellNet enables cross-study analysis of cell fate engineering protocols. Stem Cell Reports. 2023 Aug 8;18(8):1721-1742. doi: 10.1016/j.stemcr.2023.06.008. Epub 2023 Jul 20. PMID: 37478860; PMCID: PMC10444577.
• eSNPKaryotyping: Weissbein U, Schachter M, Egli D, Benvenisty N. Analysis of chromosomal aberrations and recombination by allelic bias in RNA-Seq. Nat Commun. 2016 Jul 7;7:12144. doi: 10.1038/ncomms12144. PMID: 27385103; PMCID: PMC4941052.
• GATK, Docker, WDL: Van der Auwera GA & O'Connor BD. (2020). Genomics in the Cloud: Using Docker, GATK, and WDL in Terra (1st Edition). O'Reilly Media.
• fastqc: Li, B., Dewey, C.N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011). https://doi.org/10.1186/1471-2105-12-323
• RSEM: Li, B., Dewey, C.N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011). https://doi.org/10.1186/1471-2105-12-323
• STAR Alignment: Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013 Jan 1;29(1):15-21. doi: 10.1093/bioinformatics/bts635. Epub 2012 Oct 25. PMID: 23104886; PMCID: PMC3530905.
• parallel: O. Tange (2018): GNU Parallel 2018, March 2018, https://doi.org/10.5281/zenodo.1146014.
• Bowtie: Langmead, B., Salzberg, S. Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357–359 (2012). https://doi.org/10.1038/nmeth.1923