Epigenetic silencing and genome dynamics determine the fate of giant virus endogenizations in Acanthamoeba

This is a repository of data files and code used and produced in the bioinformatic analyses for the publication "Epigenetic silencing and genome dynamics determine the fate of giant virus endogenizations in Acanthamoeba." Data Genome assembly and gene predictions The genome assemblies and gene models that form the backbone of this publication were generated as part of work for the Acanthamoeba genome paper by Matthey-Doret, Colp et al. (2022). Because the two projects ran partly in parallel, the assemblies and gene models used here represent different stages of the work by Matthey-Doret, Colp et al. (2022), and do not perfectly match the records publically provided for this publication. As such, we provide here the version of the files used for this study. Archive: Assemblies_and_genes_Colp_Matthey-Doret.tar.gz Protein sequences and function This directory includes fastas for predicted protein (based on the gene models discussed above) and intergenic ORFs (identifiable from the "extraction" in their ID). Interproscan files for both Neff and C3 are also included. Archive: Proteins_and_function.tar.gz Viral genes, viral regions, divergent regions This repository holds the relevant coordinates and information on viral genes. The "*_all_genes_info_virus_conservation.bed" files have the coordinates for all genes and specifies whether they are viral or non viral, and whether they are conserved, degraded or unconserved (for the purpose of all analyses, degraded genes are treated as conserved). The "*_mummer_concatenated_divergent_regions.gff" files have the coordinates of all divergent regions in Neff and C3, and specifies which ones contain all least one full viral gene. The "*_viral_regions_curated.gff" files show the coordinates for manually curated viral regions based on the divergent regions. Archive: Viral_annotations.tar.gz Genome alignment information This repository holds the output files for the mummer4 Neff-C3 genome alignments. Archive: mummer4_genome_alignments.tar.gz Methylation Methylation frequency at every CpG site in Neff and C3. Calls are based on raw reads provided by Colp and Matthey-Doret. Due to the large size of raw nanopore fast5 files, the raw nanopore files are only available on demand. Archive: Methylation.tar.gz Custom viral regions python script input The customized .gff files indicating whether a gene is viral are in the format necessary to be used with the script mummer_concat_for_publication_inputs.py. Archive: input_for_viral_regions_script.tar.gz BLAST hits and taxonomy The initial output file for the Neff and C3 diamond blast searches (diamond_*_v_NCBI.tsv), modified to include the taxonomy of subject proteins, as well as the filtered version including only viral candidates (diamond_*_viral_candidates.tsv). Files breaking down the taxonomy of viral proteins are also included. Archive: BLAST.tar.gz Viral hallmark genes Output files for the ViralRecall hmm viral mhallmark gene search and selected list of viral hallmark genes found in viral regions, which were used for figure 2. Archive: ViralRecall_v2.tar.gz Expression Transcriptome fasta for Neff and C3, based on the exon sequences of the Colp and Matthey-Doret gene models and the nucleotide sequence of viral intergenic ORFs. Kallisto output files are also included. These files are modified to include information about which genes are viral or not and conserved or not. Archive: Expression.tar.gz Mobile elements Repeatmasker .out files, with simple and low complexity repeats filtered out, and modified to be easier to process in R. Archive: Mobile_elements.tar.gz Trees Alignment and tree files for all phylogenetic trees featured in the study. Archive: Trees.tar.gz Macrosynteny .bed coordinates of all protein-coding sequences used in synteny analysis and orthologue table used for the same analysis. Archive: Macrosynteny.tar.gz Circos plot A manually curated set of viral DNA sequences was used to generate the circos alignment in fig. 3. The script was run with circos_candidates_v3.fa, although a version with more instructive fasta headers can be found at circos_viral_regions.fa. The initial alignment and configuration files were generated using Sibelia, and the image was rendred based on the configuration files in the directory based on Circos. Several annotation files were created to show the location of genomic features, and a custom python script is included in the directory which re-colours links based on the genomic features they overlap. This is the only script included in the data rather than code section, as its extremely specific use case doesn't make sense outside of the exact directory it was created for. Note that some configuration files were further manually edited to achieve the desired aesthetic, and that the final image was further edited in Affinity Designer. Abandon all hope all yee who enter here. Archive: custom_circos.tar.gz Distances from chromosome ends Files showing the distance from the nearest chromosome ends of various genomic features in Neff and C3. The distance is showed in the second to last column of each .tsv file. Distances of Neff genomic features were corrected to account for the "true" start of chromosomes, determined based on a manual assessment of misassemblies using long-reads mapping to the genome and Hi-C data; where the extremity of a scaffold was deemed to be a misassembly, distance was calculated based on the start or end provided in 12_Ac_Neff_rcor_polish_corrected_ends.tsv. Sequences overlapping with suspicious regions (the coordinates for which are found in Neff_Misassembly_filter_list.tsv) were filtered out of the analysis. While the same filtering step was not done for C3, "null" filtering files are provided for C3 to run the scripts. Archive: distances.tar.gz Code All python scripts are found in: custom_python.tar.gz All R scripts are found in: R_scripts.tar.gz All bash scripts are found in: program_submission_scripts.tar.gz Viral sequence detection Viral proteins were detected using a diamond BLASTp search of all predicted proteins and intergenic ORFs in Neff and C3 against a local copy of the NCBI database. qsub diamond_blast.sh Taxonomy was assigned to hits using acc2tax. Hits to other Acanthamoeba proteins, as well as entries with taxonomy listed as 'NO TAXA FOUND', 'other entries' and 'Unknown' were filtered out using grep -v ',Acanthamoeba,' $ | grep -v 'NO TAXA FOUND' | grep -v 'other entries' | grep -v 'Unknown'. Proteins were retained as viral if either the top hit after filtering was viral, or if at least 5/10 top blast hits were viral. Viral hallmark genes were detected using ViralRecall: qsub viral_recall_hallmark_genes.sh Conservation and delineation of viral regions The Neff and C3 genomes were aligned with mummer4, and the show-diff function was used to identify differences in the assemblies: qsub mummer4_genome_alignment_final.sh Overlapping differences in the alignments were concatenated together to form 'divergent regions', using a custom script which accounts for translocations. Viral regions were defined as divergent regions which had at least one fully overlapping gene. Unconserved genes were defined as genes which fully overlapped with show-diff features corresponding to gaps in the alignment (GAP, BRK, JMP, INV, SEQ). python mummer_concat_for_publication_inputs.py Expression To assess global expression levels, HISAT2 was used to align reads to the reference: qsub hisat2_C3_final.sh The specific transcripts per million (TPM) values for all predicted genes and viral intergenic ORFs were calculated using kallisto: qsub kallisto_final.sh for Neff. qsub kallisto_single_final.sh for c3. A t-test comparing the TPM values of conserved and unconserved viral and non-viral genes was performed: Rscript expression_violin.R Methylation analysis Methylation calls using nanopolish require aligning nanopore reads to a genome reference. This was done with minimap2: qsub minimap2_final.sh The resulting .bam file is then used as an input file for the nanopolish script: qsub call_methylation_nanopolish_final.sh A t-test comparing the methylation levels of CpG sites across different genomic contexts was then performed: Rscript genomic_context_methylation_violin.R Mobile elements Mobile elements were predicted using repeatmasker: qsub repeat_masker.sh Low_complexity and Simple_repeat predictions were removed from the repeatmasker .out file using grep -v $file.out. The file was then made R friendly by replacing sequences of multiple spaces with tabs and removing leading characters: sed 's/ */ /g' file.out | cut -c2- | sed 's/ / /g' A chi-square test was performed to compare the statistical significance of differences in mobile element composition of different genomic contexts: Rscript chi_square_repeat_masker.R Functional annotation The function of Neff proteins was predicted using interposcan: qsub interproscan.sh Phylogenetic trees Phylogenetic trees were generated using IQ-tree: qsub iqtree.sh Distance from chromosome ends The scripts used to calculate the distance of genomic features from the nearest chromosome end in the first 32 scaffods of Neff and C3 (which were deemed large enough to be "chromosome scale" for the purpose of our analsis) can be run with: python get_distance_curated_C3_for_publication.py python get_distance_curated_Neff_for_publication.py Macrosynteny The macrosynteny analysis can be run with the following R script: Rscript Acas_macrosynteny.R Circos plot The alignment underlying the circos plot in fig. 3 was generated using the following Sibelia command: Sibelia -s loose -m 500 circos_candidates_v3.fa And then calling the circos program from the generated circos directory, which contained painstakingly made annotations files with the location of overlaping genomic features, as well as modified configuration files to obtain the correct aesthetic, and a custom python script to recolor connections based on the genomic features they overlap (found in custom_circos.tar.gz).本仓库收录了为发表论文《表观遗传沉默与基因组动态调控决定棘阿米巴中巨型病毒内源化的命运》所开展的生物信息学分析中使用及生成的数据文件与代码。 ## 数据 ### 基因组组装与基因预测 本论文核心依托的基因组组装结果与基因模型，由Matthey-Doret、Colp等人（2022）发表的棘阿米巴基因组研究工作中生成。由于两项研究部分并行开展，本研究使用的组装结果与基因模型对应Matthey-Doret、Colp等人（2022）工作的不同阶段，与该论文公开的记录不完全一致。因此，本仓库提供本研究实际使用的文件版本。 存档文件：Assemblies_and_genes_Colp_Matthey-Doret.tar.gz ### 蛋白质序列与功能注释 该目录包含基于上述基因模型预测得到的蛋白质序列，以及基因间区开放阅读框（intergenic ORFs，可通过其ID中的“extraction”字段识别）的FASTA格式文件。同时收录了Neff与C3两个样本的InterProScan分析结果文件。 存档文件：Proteins_and_function.tar.gz ### 病毒基因、病毒区域与差异区域 本仓库收录了病毒基因的相关坐标与注释信息。其中，`*_all_genes_info_virus_conservation.bed`文件包含所有基因的坐标，并标注其为病毒基因或非病毒基因，以及是否保守、降解或未保守（本研究所有分析中，降解基因均视为保守基因）。`*_mummer_concatenated_divergent_regions.gff`文件包含Neff与C3基因组中所有差异区域的坐标，并标注其中包含至少一个完整病毒基因的区域。`*_viral_regions_curated.gff`文件则展示了基于差异区域手动整理的病毒区域坐标。 存档文件：Viral_annotations.tar.gz ### 基因组比对信息 本仓库收录了mummer4比对Neff与C3基因组得到的输出文件。 存档文件：mummer4_genome_alignments.tar.gz ### 甲基化分析 收录Neff与C3样本中所有CpG位点的甲基化频率数据。甲基化位点识别基于Colp与Matthey-Doret提供的原始测序reads。由于纳米孔（nanopore）测序原始fast5格式文件体积较大，仅可按需获取。 存档文件：Methylation.tar.gz ### 病毒区域自定义Python脚本输入文件 标注基因是否为病毒基因的自定义GFF格式文件，可直接用于脚本`mummer_concat_for_publication_inputs.py`。 存档文件：input_for_viral_regions_script.tar.gz ### BLAST比对结果与分类学注释 收录Neff与C3样本的Diamond BLASTp比对初始输出文件（diamond_*_v_NCBI.tsv），该文件已补充比对蛋白的分类学注释；同时收录仅保留病毒候选序列的过滤版本文件（diamond_*_viral_candidates.tsv）。此外还包含病毒蛋白分类学信息拆分文件。 存档文件：BLAST.tar.gz ### 病毒标志性基因 收录ViralRecall软件的隐马尔可夫模型（Hidden Markov Model, HMM）病毒标志性基因搜索输出文件，以及在病毒区域中鉴定得到的标志性病毒基因列表（用于绘制图2）。 存档文件：ViralRecall_v2.tar.gz ### 转录组表达分析 收录基于Colp与Matthey-Doret基因模型的外显子序列，以及病毒基因间区开放阅读框的核苷酸序列生成的Neff与C3转录组FASTA文件。同时包含Kallisto转录本定量输出文件，所有文件均已补充基因是否为病毒基因、是否保守的注释信息。 存档文件：Expression.tar.gz ### 可移动元件分析 收录经过过滤的RepeatMasker重复序列注释软件输出文件（.out格式），已移除简单重复与低复杂度重复序列，并进行格式调整以适配R语言处理。 存档文件：Mobile_elements.tar.gz ### 系统发育树分析 收录本研究中所有系统发育树的比对文件与树文件。 存档文件：Trees.tar.gz ### 宏观共线性分析 收录共线性分析中使用的所有蛋白质编码序列的BED格式坐标文件，以及该分析所用的同源基因表。 存档文件：Macrosynteny.tar.gz ### Circos可视化图谱 本研究使用手动整理的病毒DNA序列集合生成图3的Circos可视化图谱。脚本基于`circos_candidates_v3.fa`运行，而带有更清晰FASTA序列头的版本可在`circos_viral_regions.fa`中获取。初始比对文件与配置文件通过Sibelia共线性比对软件生成，后续基于Circos软件与本目录中的配置文件绘制图像。本目录包含多个标注基因组特征位置的注释文件，以及一个自定义Python脚本，可根据比对连线重叠的基因组特征调整连线颜色。本仓库仅在此数据部分收录该脚本，因其使用场景高度特定，仅可在其生成的专属目录中正常运行。请注意，部分配置文件已进行手动编辑以达到预期的可视化效果，最终图像亦在Affinity Designer中进行了后期调整。望诸君慎入。 存档文件：custom_circos.tar.gz ### 染色体末端距离分析 收录Neff与C3样本中各类基因组特征距最近染色体末端的距离文件，距离信息位于每个TSV文件的倒数第二列。Neff样本的基因组特征距离已进行校正，以适配染色体的“真实”起始位置：该校正基于长读长序列比对基因组与Hi-C数据手动评估的组装错误结果；若某支架的末端被判定为组装错误，则距离计算将基于`12_Ac_Neff_rcor_polish_corrected_ends.tsv`中提供的起始或终止位置。与可疑区域（坐标见`Neff_Misassembly_filter_list.tsv`）重叠的序列已从分析中过滤。C3样本未执行相同的过滤步骤，但提供了空过滤文件以支持脚本运行。 存档文件：distances.tar.gz ## 代码 所有Python脚本收录于`custom_python.tar.gz`；所有R脚本收录于`R_scripts.tar.gz`；所有Bash脚本收录于`program_submission_scripts.tar.gz`。 ### 病毒序列鉴定 病毒蛋白鉴定通过对Neff与C3样本中所有预测蛋白及基因间区开放阅读框进行Diamond BLASTp比对，比对数据库为本地NCBI蛋白库。比对任务通过`qsub diamond_blast.sh`提交。使用`acc2tax`工具为比对结果分配分类学注释。通过以下命令过滤掉比对到其他棘阿米巴蛋白的结果，以及分类学标注为“NO TAXA FOUND”“other entries”与“Unknown”的条目：`grep -v ',Acanthamoeba,' $ | grep -v 'NO TAXA FOUND' | grep -v 'other entries' | grep -v 'Unknown'`。若过滤后的最优比对结果为病毒蛋白，或前10个最优比对结果中有至少5个为病毒蛋白，则该蛋白被判定为病毒蛋白。病毒标志性基因通过ViralRecall软件鉴定：`qsub viral_recall_hallmark_genes.sh` ### 病毒区域的保守性界定 使用mummer4软件比对Neff与C3基因组，并通过`show-diff`功能识别组装差异：`qsub mummer4_genome_alignment_final.sh`。利用考虑了易位事件的自定义脚本，将比对中重叠的差异区域拼接为“差异区域”。将包含至少一个完整重叠基因的差异区域定义为病毒区域。将与比对间隙（GAP、BRK、JMP、INV、SEQ）对应的`show-diff`特征完全重叠的基因定义为非保守基因。运行脚本：`python mummer_concat_for_publication_inputs.py` ### 转录组表达分析 为评估全局表达水平，使用HISAT2测序比对软件将测序reads比对至参考基因组：`qsub hisat2_C3_final.sh`。使用Kallisto转录本定量软件计算所有预测基因与病毒基因间区开放阅读框的转录本每百万reads计数（Transcripts Per Million, TPM）：Neff样本对应任务为`qsub kallisto_final.sh`，C3样本对应任务为`qsub kallisto_single_final.sh`。通过`t`检验比较保守与非保守的病毒基因及非病毒基因的TPM值：`Rscript expression_violin.R` ### 甲基化分析 使用nanopolish甲基化分析软件进行甲基化位点识别需先将纳米孔测序reads比对至参考基因组，该步骤通过minimap2比对软件完成：`qsub minimap2_final.sh`。生成的BAM格式文件将作为输入文件用于nanopolish脚本：`qsub call_methylation_nanopolish_final.sh`。通过`t`检验比较不同基因组背景下CpG位点的甲基化水平：`Rscript genomic_context_methylation_violin.R` ### 可移动元件分析 可移动元件预测通过RepeatMasker软件完成：`qsub repeat_masker.sh`。通过以下命令从RepeatMasker的.out输出文件中移除低复杂度重复与简单重复注释：`grep -v $file.out`。随后通过以下命令调整文件格式以适配R语言处理：将多个空格替换为制表符并移除首字符：`sed 's/ */ /g' file.out | cut -c2- | sed 's/ / /g'`。通过卡方检验比较不同基因组背景下可移动元件组成的差异显著性：`Rscript chi_square_repeat_masker.R` ### 功能注释 Neff样本的蛋白功能通过InterProScan软件预测：`qsub interproscan.sh` ### 系统发育树分析 系统发育树通过IQ-tree系统发育树构建软件生成：`qsub iqtree.sh` ### 染色体末端距离分析 用于计算Neff与C3前32个scaffolds（经评估其长度足够达到“染色体级别”）中各类基因组特征距最近染色体末端距离的脚本，可通过以下命令运行：`python get_distance_curated_C3_for_publication.py`与`python get_distance_curated_Neff_for_publication.py` ### 宏观共线性分析 宏观共线性分析可通过以下R脚本运行：`Rscript Acas_macrosynteny.R` ### Circos可视化图谱 图3中Circos图谱对应的比对结果通过以下Sibelia命令生成：`Sibelia -s loose -m 500 circos_candidates_v3.fa`。随后在生成的Circos目录中调用Circos软件，该目录包含精心制作的基因组特征重叠位置注释文件、用于调整可视化效果的自定义配置文件，以及一个可根据重叠基因组特征调整连线颜色的自定义Python脚本（收录于`custom_circos.tar.gz`）。