Processing Data_Nina_RNAseq_2026_dataset_2

Leave a reply

一个重要的区别是，在本次实验中，我们使用了不同的VZV毒株（EMC-1）来感染类器官。 我附上了一份包含该序列的GenBank文件。关于这个参考序列，有几点需要注意：

  • 它是完整的BAC构建体序列（135 kb），而不是单纯的病毒基因组。它包含了BAC骨架（mini-F、cat/氯霉素抗性基因、loxP等），这些在受感染的类器官中不应该被转录，此外还包含一个EGFP表达盒。请注意，GFP报告基因在受感染的细胞中是有表达的，因此数据中会存在GFP的测序reads——请在参考序列中保留GFP序列，以便这些reads能够正确比对，而不是丢失或发生错误比对。

  • 序列本身是EMC-1的序列，但其注释是基于相似性从Dumas参考序列（NC_001348）转移过来的，因此ORF（开放阅读框）的坐标可能并不完全准确。重复区域（TRL/IRL、TRS/IRS以及R3/R4重复序列）的注释可能不太可靠，并且由于IRS/TRS的结构，ORF62/63/64会出现两次。

  • 该文件中包含几个已注释的变异位点。这些大多是预期的毒株差异（EMC-1与Dumas相比），但也有少数可能是组装伪影（assembly artifacts），所以我不能完全确定每一个碱基对于我们实际的病毒来说都是绝对正确的，因为这个序列文件我是从另一个实验室拿到的。

Based on the pipeline established in your previous analysis and Nina’s specific requirements for the new EMC-1 VZV strain, here is the fully adapted workflow.

I have specifically addressed her concerns regarding the EGFP cassette, the duplicated ORF62/63/64 regions, and the Salmon requirements for exon features.

Step 1: Prepare samplesheet.csv

Create a samplesheet.csv in your working directory (/mnt/md1/DATA/Data_Nina_RNAseq_2026_batch_2/). The sample names are kept consistent with your previous batch to ensure the R scripts will work seamlessly.

sample,fastq_1,fastq_2,strandedness
control_r1,./20261020_AV243904_0124_B_CORALL_Nina/01_negative_42_91/01_negative_42_91_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/01_negative_42_91/01_negative_42_91_R2.fastq.gz,auto
control_r2,./20261020_AV243904_0124_B_CORALL_Nina/02_negative_42_93/02_negative_42_93_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/02_negative_42_93/02_negative_42_93_R2.fastq.gz,auto
control_r3,./20261020_AV243904_0124_B_CORALL_Nina/03_negative_42_94/03_negative_42_94_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/03_negative_42_94/03_negative_42_94_R2.fastq.gz,auto
VZV.d10_r1,./20261020_AV243904_0124_B_CORALL_Nina/04_VZV_d10_42_99/04_VZV_d10_42_99_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/04_VZV_d10_42_99/04_VZV_d10_42_99_R2.fastq.gz,auto
VZV.d10_r2,./20261020_AV243904_0124_B_CORALL_Nina/05_VZV_d10_42_100/05_VZV_d10_42_100_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/05_VZV_d10_42_100/05_VZV_d10_42_100_R2.fastq.gz,auto
VZV.d10_r3,./20261020_AV243904_0124_B_CORALL_Nina/06_VZV_d10_42_101/06_VZV_d10_42_101_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/06_VZV_d10_42_101/06_VZV_d10_42_101_R2.fastq.gz,auto
VZV.d15_r1,./20261020_AV243904_0124_B_CORALL_Nina/07_VZV_d15_42_105/07_VZV_d15_42_105_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/07_VZV_d15_42_105/07_VZV_d15_42_105_R2.fastq.gz,auto
VZV.d15_r2,./20261020_AV243904_0124_B_CORALL_Nina/08_VZV_d15_42_106/08_VZV_d15_42_106_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/08_VZV_d15_42_106/08_VZV_d15_42_106_R2.fastq.gz,auto
VZV.d15_r3,./20261020_AV243904_0124_B_CORALL_Nina/09_VZV_d15_42_113/09_VZV_d15_42_113_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/09_VZV_d15_42_113/09_VZV_d15_42_113_R2.fastq.gz,auto
VZV.d20_r1,./20261020_AV243904_0124_B_CORALL_Nina/10_VZV_d20_42_118/10_VZV_d20_42_118_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/10_VZV_d20_42_118/10_VZV_d20_42_118_R2.fastq.gz,auto
VZV.d20_r2,./20261020_AV243904_0124_B_CORALL_Nina/11_VZV_d20_42_122/11_VZV_d20_42_122_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/11_VZV_d20_42_122/11_VZV_d20_42_122_R2.fastq.gz,auto
VZV.d20_r3,./20261020_AV243904_0124_B_CORALL_Nina/12_VZV_d20_42_126/12_VZV_d20_42_126_R1.fastq.gz,./20261020_AV243904_0124_B_CORALL_Nina/12_VZV_d20_42_126/12_VZV_d20_42_126_R2.fastq.gz,auto

Step 2: Robust GenBank to FASTA/GTF Conversion

Save the following Python script as gb_to_fasta_gtf_emc1.py. It specifically addresses Nina’s notes:

  1. Preserves EGFP: It processes all features, so the EGFP cassette and BAC backbone will be naturally included.
  2. Handles Duplicated ORFs: It automatically appends _copy2, _copy3, etc., to duplicated transcript IDs (like ORF62/63/64 in the IRS/TRS regions) to prevent Salmon from crashing on duplicate IDs.
  3. Synthesizes exon features: As noted in your previous DEBUG_2, Salmon strictly requires exon features. If the GenBank lacks them, this script safely synthesizes them from the CDS or gene boundaries.
#!/usr/bin/env python3
import sys
from Bio import SeqIO

if len(sys.argv) != 4:
    print("Usage: python3 gb_to_fasta_gtf_emc1.py <input.gbk> <output.fasta> <output.gtf>")
    sys.exit(1)

gb_file = sys.argv[1]
fasta_file = sys.argv[2]
gtf_file = sys.argv[3]

try:
    record = SeqIO.read(gb_file, "genbank")
except Exception as e:
    print(f"Error reading GenBank file: {e}")
    print("Note: If this is a date format error (e.g. DD-MM-YYYY instead of DD-MMM-YYYY), fix it using sed before running this script.")
    sys.exit(1)

seq_id = record.id
if not seq_id or seq_id == ".":
    seq_id = "EMC1_BAC"

# 1. Write FASTA
with open(fasta_file, "w") as f:
    SeqIO.write(record, f, "fasta")
print(f"Successfully wrote {fasta_file}")

# 2. Write GTF
with open(gtf_file, "w") as f:
    has_exon = any(feat.type == "exon" for feat in record.features)
    seen_transcripts = {} # To handle duplicated ORFs (e.g., ORF62/63/64)

    for feature in record.features:
        if feature.type in ["gene", "mRNA", "transcript", "CDS", "exon"]:
            # Prioritize 'gene', then 'locus_tag'
            gene_id_candidates = feature.qualifiers.get("gene", []) + feature.qualifiers.get("locus_tag", [])
            gene_id = str(gene_id_candidates[0]) if gene_id_candidates else "unknown_gene"
            gene_id = gene_id.replace(" ", "_").replace('"', '')

            transcript_id_candidates = feature.qualifiers.get("transcript_id", [])
            transcript_id = str(transcript_id_candidates[0]) if transcript_id_candidates else f"{gene_id}_RNA"
            transcript_id = transcript_id.replace(" ", "_").replace('"', '')

            # Handle duplicates (e.g., ORF62 appearing twice in IRS/TRS)
            if transcript_id in seen_transcripts:
                seen_transcripts[transcript_id] += 1
                transcript_id_unique = f"{transcript_id}_copy{seen_transcripts[transcript_id]}"
            else:
                seen_transcripts[transcript_id] = 1
                transcript_id_unique = transcript_id

            start = feature.location.start + 1
            end = feature.location.end
            strand = "+" if feature.location.strand == 1 else "-"

            # Build attributes
            attrs = [f'gene_id "{gene_id}"', f'transcript_id "{transcript_id_unique}"']
            if "product" in feature.qualifiers:
                prod = str(feature.qualifiers["product"][0]).replace('"', "'")
                attrs.append(f'gene_name "{prod}"')
            attr_str = "; ".join(attrs) + ";"

            # Write the main feature
            f.write(f"{seq_id}\tGenBank\t{feature.type}\t{start}\t{end}\t.\t{strand}\t.\t{attr_str}\n")

            # Salmon requires 'exon' features. If missing, synthesize from CDS or gene
            if not has_exon and feature.type in ["gene", "CDS", "mRNA", "transcript"]:
                exon_attr_str = f'gene_id "{gene_id}"; transcript_id "{transcript_id_unique}"; exon_id "{transcript_id_unique}_exon1";'
                f.write(f"{seq_id}\tGenBank\texon\t{start}\t{end}\t.\t{strand}\t.\t{exon_attr_str}\n")

print(f"Successfully wrote {gtf_file}")

Run the conversion:

cd /mnt/md1/DATA/Data_Nina_RNAseq_2026_batch_2
# Optional: Fix date format if Biopython complains about DD-MM-YYYY
# sed -i 's/13-06-2025/13-JUN-2025/' EMC1BACreGFPT2AORF4_1.gbk 

python3 gb_to_fasta_gtf_emc1.py EMC1BACreGFPT2AORF4_1.gbk EMC1_final.fasta EMC1_final.gtf

Step 3: Nextflow Execution Commands

Run the pipeline in two parts (Host and Virus). The parameters remain identical to your previous run to maintain consistency.

A. Host Analysis (GRCh38)

nextflow run /home/jhuang/Tools/nf-core-rnaseq-3.12.0/main.nf -profile docker \
--input samplesheet.csv \
--outdir results_GRCh38 \
--genome GRCh38 \
--with_umi \
--umitools_extract_method regex \
--umitools_bc_pattern '^(?P
<umi_1>.{12}).*' \
--trimmer fastp \
--extra_fastp_args "--trim_front2 10" \
-resume \
--max_cpus 54 \
--max_memory 120.GB \
--max_time 2400.h \
--save_align_intermeds \
--save_unaligned \
--save_reference \
--aligner star_salmon \
--pseudo_aligner salmon \
--umitools_grouping_method unique

B. Virus Analysis (EMC-1) Note: --min_mapped_reads 0 is retained to prevent the pipeline from crashing on the negative control samples, which will have near-zero viral reads.

nextflow run /home/jhuang/Tools/nf-core-rnaseq-3.12.0/main.nf -profile docker \
--input samplesheet.csv \
--outdir results_EMC1 \
--fasta EMC1_final.fasta \
--gtf EMC1_final.gtf \
--with_umi \
--umitools_extract_method regex \
--umitools_bc_pattern '^(?P
<umi_1>.{12}).*' \
--trimmer fastp \
--extra_fastp_args "--trim_front2 10" \
-resume \
--max_cpus 54 \
--max_memory 120.GB \
--max_time 2400.h \
--save_align_intermeds \
--save_unaligned \
--save_reference \
--aligner star_salmon \
--pseudo_aligner salmon \
--gtf_extra_attributes gene_id \
--gtf_group_features transcript_id \
--featurecounts_group_type gene_id \
--featurecounts_feature_type transcript \
--skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --skip_multiqc \
--umitools_grouping_method unique \
--min_mapped_reads 0

(Don’t forget to fix the collated_versions.yml matplotlib bug afterwards if it appears, as noted in your previous README).

Step 4: R Script Adaptation

In your R script for merging and DEG analysis, you only need to update the setwd() paths to point to the new _batch_2 directories. The sample names and logic remain identical.

# ------ For human genome ------
setwd("~/DATA/Data_Nina_RNAseq_2026_batch_2/results_GRCh38/star_salmon")
# ... (rest of the human genome code remains the same) ...

# ------ For virus genome ------
setwd("~/DATA/Data_Nina_RNAseq_2026_batch_2/results_EMC1/star_salmon")
# ... (rest of the virus genome code remains the same) ...

# ------ Merge the raw counts of human and microbe ------
# Update the paths in the bash commands or R read.csv to point to the new batch_2 folders
# REMINDER: Manually set all counts for all virus transcripts in the three control-conditions ("control_r1","control_r2","control_r3") as 0 using Excel before merging!

Step 5: Draft Reply Email to Nina

Here is a professional and reassuring draft you can send to Nina. It explicitly confirms that her technical concerns have been addressed.

Subject: Re: Bulk RNA-seq VZV (EMC-1) infected skin organoids / Analysis Plan

Thank you for reaching out and providing the details for the new VZV RNA-seq dataset. I have reviewed the parameters and the specific notes regarding the EMC-1 strain, and I am ready to proceed with the analysis.

To ensure the highest accuracy, I have adapted our established pipeline to specifically address the points you mentioned:

  1. EGFP & BAC Backbone: The reference will include the full BAC construct sequence. The EGFP cassette will be retained in the annotation so that GFP reads map correctly and are not lost or mismapped.
  2. Duplicated ORFs (IRS/TRS regions): I have implemented a custom script to process the GenBank file that automatically detects and resolves duplicated transcript IDs (e.g., ORF62/63/64 appearing twice) by appending unique identifiers. This prevents the quantification software (Salmon) from crashing on duplicate entries.
  3. Annotation Gaps: The script also synthetically generates the mandatory exon features required by the pipeline if they are missing from the transferred Dumas annotation, ensuring robust mapping.
  4. Negative Controls: As with the previous batch, the pipeline is configured to gracefully handle the negative control samples, which will naturally have near-zero viral reads, without aborting the workflow.

The sequencing parameters (UMI, 2x150bp, Lexogen CORALL v2) match the previous run perfectly, so the same UMI extraction and trimming logic will be applied.

I will begin processing the data from the provided LIV link shortly. I expect to have the initial QC, PCA plots, and differential expression results (Host + Viral) ready for your review within the next few days.

Please let me know if you have any additional preferences or if there is anything else I should keep in mind.

Leave a Reply

Your email address will not be published. Required fields are marked *