本期教程
为方便大家的学习，我们提供此文的原文及译文（机器翻译）和分析代码，在后台回复关键词**：20240709，即可获得。请大家看清楚回复关键词，每天都有很多人回复错误关键词，我这边没时间和精力一一回复。**
写在前面
在2023年4月6日，在Nature Genetics中发表了番茄的超级泛基因的文章，Super-pangenome analyses highlight genomic diversity and structural variation across wild and cultivated tomato species,该文章使用的9个野生种和2个栽培品种，阐明了茄属番茄组（Solanum section Lycopersion）的基因组演化历史，构建了首个番茄超级泛基因组/图基因组。
时隔一年后，2024年7月8日，北京大学现代农业研究院在Nature Genetics发表西瓜的超级泛基因组文章Telomere-to-telomere Citrullus Super-pangenome Provides Direction for Watermelon Breeding，完成了27个不同西瓜基因型的端粒到端粒(T2T)组装。
该研究内容获得众多**「育种科学家」**的点评，是一个非常完美的工作。
摘要
为了破译瓜属西瓜的遗传多样性，我们获得七个西瓜品种，共27个不同基因型的端粒到端粒(T2T)组装。这个T2T超级全基因组通过添加399.2Mb 和11,225个基因扩展了先前发表的参考基因组T2T-G42。比较分析揭示了西瓜的基因变异和结构变异(SVs) ，揭示了西瓜的进化和驯化过程，这些过程提高了西瓜的苦味和含糖量等属性，同时降低了西瓜的抗病性。西瓜和西瓜的多重抗病基因座被成功地引入栽培西瓜。在西瓜中发现的SVs不仅遗传自科尔多芬属，而且也遗传自粘籽隐杆线虫，这表明在栽培西瓜的谱系中还有其他科尔多芬属以外的祖先。我们的研究大大提高了对西瓜基因组多样性的理解，为所有西瓜品种提供了全面的参考基因组。这一进步有助于利用西瓜的野生亲缘基因治疗进行西瓜的开发和利用。
原文内容
以下内容来自公众号“BioArt植物**”在2024年7月8日发表的推文，原文链接：**https://mp.weixin.qq.com/s/Rl4hA5o3MNFLdMzM8xYaBg
泛基因组是一个物种中所有个体基因组信息的总和，构建泛基因组可以有效解决单一参考基因组带来的信息缺失和分析偏差。而超级泛基因组则代表一个属内所有物种的基因组信息，尤其蕴含了野生种中丰富的基因组变异，是对泛基因组的进一步扩展，在远缘杂交和基因发掘等方向具有重要应用前景。
西瓜是世界重要的园艺经济作物之一，富含抗氧化剂番茄红素和增强血液循环的瓜氨酸，口味甘甜，享有夏季"水果之王"的美誉，深受全球消费者喜爱。西瓜全球年产量接近1亿吨，我国作为种瓜和吃瓜第一大国，生产和消费了全球总量的60%以上，西瓜产业在乡村振兴和农民致富中占有十分重要的作用。我国西瓜产业的迅速发展离不开优良的品种 ，依靠一代代育种家的不懈努力，我国成为世界上西瓜种质创新和品种选育最为活跃的国家，在品种种源方面实现了完全自主可控。随着气候变化的加剧和病虫危害的日益严重，西瓜生产面临着严重的挑战。在现有优良品种的基础上，“植入”缺少的抗病耐逆优良基因，是西瓜种质创新的关键问题，对维持我国西瓜产业高质量绿色发展十分重要。栽培西瓜在驯化和品种改良过程中因为对品质和产量的追求，导致遗传多样性非常狭窄，抗病耐逆基因大量丢失，而在西瓜的祖先即野生西瓜中存在广泛的遗传与表型多样性，具有丰富的抗病耐逆等基因资源，是西瓜遗传改良和种质创新的宝库。
该研究绘制了西瓜属全部7个种的28份代表性材料的端粒到端粒（T2T）高质量基因组图谱，成功构建首个属级T2T水平超级泛基因组。得到总计768.5Mb，32,513个基因家族的西瓜属泛基因组，是单个西瓜基因组的1.5倍，增加了11,225个栽培西瓜中没有发现的基因。基于T2T高质量基因组图谱，该研究首次全面比较了西瓜属的着丝粒序列，其丰富的变异和特有的进化关系影响了不同种之间的杂交配对。因此，该超级泛基因组极大地扩展了西瓜遗传改良的基因池（图1）。
图1：西瓜属7个种的28份代表性材料的多样性与泛基因组图谱
利用基因组序列，该研究拓展了西瓜属的分类系统，核实了西瓜起源于非洲的理论依据，并发现栽培西瓜除之前报道的cordophanus亚种外可能还存在其他祖先。另外，该研究还鉴定到在西瓜属内发生的三次重大染色体重排事件和两个超长片段倒位，这些染色体重排显著影响了西瓜的抗性、品质相关基因以及三维基因组结构，并在栽培种中得以保留。同时，这些发现也为回交育种过程中避免连锁累赘提供了基因组基础（图2）。
图2：西瓜属进化与栽培西瓜的起源
西瓜超级泛基因组鉴定出461,987个结构变异（SV），构建了西瓜图形泛基因组。并通过SV-GWAS鉴定了驯化过程中丢失和获得的关键基因，挖掘了与葫芦素含量、含糖量、果肉着色等重要性状相关的功能基因结构变异，发现在西瓜驯化过程中，伴随着多个与含糖量增加、果肉变红相关的基因簇扩张，大量抗病功能相关的基因簇丢失。该研究为西瓜研究人员和育种家提供了最完整的基因组资源，助力深入理解西瓜基因组的复杂性和多样性，从而高效挖掘和利用野生西瓜种中的优异基因（图3）。
图3：西瓜属重要性状基因结构变异图谱
最后，该研究利用野生西瓜的基因组序列和抗病基因信息，通过种间杂交育成了抗多种病害的自交系‘PKR6’，并确定了西瓜抗枯萎病生理小种1的候选基因。该研究提供了利用野生种质创制优异育种材料的范例，从而将驯化过程中丢失的抗病基因重新有目的地“植入”到栽培种中，对加速抗病品种选育、促进西瓜产业高效发展具有深远意义（图4）。
图4：远缘杂交创制优良抗病种质PKR6基因组及枯萎病抗性表型
西瓜属端粒到端粒超级泛基因组是最高质量最全面的西瓜属基因组序列图谱和变异图谱，它揭示了西瓜属的基因组演化历史，发掘了西瓜野生种中丰富的遗传多样性，提供了利用野生种质创制优异育种材料的新范例，为其他作物超级泛基因组的构建和野生种质的利用指明了方向。
专家点评
黄三文 院士 点评
参考基因组是分子育种的重要基础。2008 年，多家科研单位和高校共同发起并推出了“国际西瓜基因组”计划。该计划预期首先实现西瓜基因组草图，之后更新为高质量基因组图谱，并进一步扩展至整个西瓜属。到2012 年，北京市农林科学院等单位联合发布了第一个东亚栽培西瓜“97103”的基因组草图，开启西瓜基因组研究时代（Nature Genetics, 2013）。之后，在多个中外团队共同努力下，于2019 年在Nature Genetics 杂志发布了改进的西瓜“97103”的参考基因组图谱，并提供了西瓜属414 份代表性材料的全基因组重测序变异图谱。该研究不仅明确了西瓜7 个种之间的进化关系，而且发掘了控制含糖量与瓤色等品质性状的关键基因，有助于开展相应性状改良的分子育种技术。然而，受到当时测序技术的限制，已报道的参考西瓜基因组序列中仍存在着大量缺口（Gaps），这些上万个缺口中含有丰富的基因和重复元件信息，因此，不完整的参考基因组制约了基因组学研究和育种应用。到2022 年，北京大学现代农业研究院终于完成了首个西瓜端粒到端粒（T2T）的基因组图谱组装，并构建相应的饱和EMS 突变体库 (Molecular Plant, 2022) 。同时这项研究也是葫芦科作物中属水平的第一个T2T 水平的参考基因组，这一里程碑式的工作不仅为西瓜提供了高分辨率的基因组图谱，也为葫芦科作物的后续遗传研究建立了新标准。
到了T2T 基因组时代，研究人员发现相比栽培西瓜的高质量基因组，野生西瓜的参考基因组仍很缺乏。栽培种西瓜的遗传变异非常狭窄，育种家已经将越来越多的西瓜野生种用于育种来提高栽培西瓜的抗病性和抗逆性，但高质量野生西瓜基因组序列的缺乏阻碍了野生资源的有效利用。为解决这一问题，西瓜基因组学研究必须包括更大的遗传多样性。2024 年这篇西瓜属T2T 水平超级泛基因组论文的发表，代表了科研界对西瓜属基因组理解和利用的重大飞跃。
首先，这项研究将对西瓜基因组理解从“种”的水平提升到“属”的水平。该研究涵盖了西瓜属全部的7 个种，其中包括6 个野生种和1 个栽培种，一共28 份材料的T2T 水平的基因组信息。“属”水平的泛基因组相比已报道的栽培种基因组增加了399.2 Mb 的基因组序列和11,225 个基因，因此极大丰富了西瓜的基因资源库。该研究鉴定了野生西瓜种大量新基因和结构变异，为基因功能研究和育种应用提供了宝贵的资源库。
其次，该研究加深了我们对西瓜进化和驯化的理解。通过比较西瓜属水平的基因组序列，进化中的染色体重排事件和着丝粒组成结构变化均是首次被准确鉴定，这些染色体水平的遗传变异揭示了当前西瓜遗传背景形成的复杂过程。其中对野生和栽培西瓜结构变异的比较，为西瓜的进化和驯化进程提供了新证据，从而得到栽培西瓜物种形成的新假说。
更重要的是，这项研究将极大推动西瓜抗病研究和育种改良进程。研究种的一个亮点是从两个野生种C. amarus（饲料西瓜）和C. mucosospermus（黏籽西瓜）中系统挖掘抗病基因，并利用杂交方式将目标基因整合到栽培西瓜基因组中。与野生种的远缘杂交不仅增强了栽培西瓜的抗病性，同时还引入了新的遗传多样性，对于西瓜长期可持续育种发展至关重要。
总之，首个T2T 水平的西瓜属超级泛基因组是令人兴奋的研究成果。这项研究不仅为西瓜野生和栽培种的基因组提供了更加全面的理解，发现了一大批新基因，为优良抗性基因的发掘和研究奠定了基础，而且对将来的育种发展具有深远影响，为更好利用野生资源改良育种指明了方向，有助于快速培育高产、优质、多抗的新品种。
此外，国际葫芦科作物大会终身成就奖获得者、美国农业部蔬菜实验室瓜类遗传资深研究员**Amnon Levi博士对该项研究亦给予了高度评价：“It is an amazing and most detailed and beautiful work ever in watermelon - 这是西瓜上至今一项不可思议、最详细和最完美的研究” 。**
北京大学现代农业研究院作物遗传与种质创新平台主要从事西瓜、甜瓜、黄瓜、南瓜和番茄等瓜菜作物种质创新和品种选育工作。在研究院农业组学大数据、生物技术等平台的大力支持下，组建四年多以来，先后发布了首个西瓜T2T基因组（Mol Plant, 2022）和T2T西瓜属超级泛基因组（Nat Genetics, 2024），利用创新的花粉EMS处理方法构建了国际上首个西瓜EMS突变体库（Mol Plant, 2022），搭建了西瓜枯萎病、蔓枯病等病害表型精准鉴定平台，挖掘了首个西瓜显性雄性不育、隐性抗白粉病（TAG, 2024）等重要农艺性状基因，基于以上基础研究全面开启了西瓜育种4.0时代。成功选育西瓜新品种7个，在第二届中国（海南）国际西甜瓜发展大会暨“海南蜜瓜”品牌推介活动会上共斩获7个奖项，包括西瓜4个组别中的3个一等奖。授权国家发明专利1项并实现了专利转让。与合作者江苏徐淮地区淮阴农业科学研究所共同获神农中华农业科技奖1项。团队研发成果引领了优质抗病西瓜品种潮流，为西瓜产业转型升级和乡村振兴提供了有力的科技支撑。
以上内容原文链接：https://mp.weixin.qq.com/s/Rl4hA5o3MNFLdMzM8xYaBg
Code availability
作者提供的超级泛基因组的组装流程，具体如下所示：
https://github.com/yilinZhang-bio/Watermelon-pangenome
01.Genome_survey
Here is provided the command to conduct a genomic survey utilizing next-generation sequencing data.
The software utilized can be acquired and employed as delineated on the specified GitHub repository.
Should you incorporate these software tools in your work, kindly acknowledge by citing this publications.
https://github.com/gmarcais/Jellyfish
https://github.com/schatzlab/genomescope
http://genomescope.org/
02.T2T_Genome_assembly_and_polishing
1. hifiasm -o ${Sample}.asm -t 40 ${Sample}.ccs.fasta
2. hifiasm -o ${Sample}.asm -t 128 --ul UL.fq.gz ${Sample}.ccs.fasta
3. lja -o ${Sample} --reads ${Sample}.ccs.fasta -t 24
4. hifiasm -o ${Sample}.asm -t 40 --write-paf --write-ec -l0 ${Sample}.ccs.fasta
rbsa.py --type nucmer --ref REF --scf SCF
filter.py --paf_fn PAF --bestn 20 --output OUTPUT
chr_paths.py --agp AGP --gfa GFA
ovlp2graph.py --overlap-file OVERLAP_FILE
graph2chr.py --reads-fasta-fn READS_FASTA_FN --paf-fn PAF_FN --sg-edges-list-fn SG_EDGES_LIST_FN --ctg-paths-fn CTG_PATHS_FN
5. nextDenovo run.cfg
6. HiC-Pro -c config-hicpro.txt -o Analysis -i HIC_read
7. python2 HiCPlotter/HiCPlotter.py -f iced/1000000/sample1_1000000_iced.matrix -o ${Sample} -r 1000000 -tri 1 -bed sample1/raw/1000000/sample1_1000000_abs.bed -n ${Sample} -wg 1 -chr Chr11
8.##get contig length
perl fastaDeal.pl -attr id:len ${Sample}.contigs.fa > ${Sample} .contigs.fa.len
##draw contig Hi-C heatmaps with 10*100000 (1-Mb) resolution
matrix2heatmap.py ${Sample}_HiC_100000_abs.bed ${Sample}_HiC_100000.matrix 10
##Run one round, when the contig assembly is quite good
perl endhic.pl ${Sample}.contigs.fa.len ${Sample}_HiC_100000_abs.bed ${Sample}_HiC_100000.matrix ${Sample}_HiC_100000_iced.matrix
ln Round_A.04.summary_and_merging_results/z.EndHiC.A.results.summary.cluster* ./
##convert cluster file to agp file
perl cluster2agp.pl Round_A.04.summary_and_merging_results/z.EndHiC.A.results.summary.cluster ${Sample}.contigs.fa.len > ${Sample}.scaffolds.agp
##get final scaffold sequence file
perl agp2fasta.pl ${Sample}.scaffolds.agp ${Sample}.contigs.fa > ${Sample}.scaffolds.fa
9. ragtag.py scaffold $REF $SCF -o ./ragtag_default -t 26
ragtag.py scaffold $REF $SCF -o ./ragtag_nucmer -t 26 --aligner 'nucmer'
10. ragtag.py patch $target $query
11.
## Screen out sequences containing telomere repeats
## For ont reads
minimap2 -t80 -x map-ont ref_genome.fasta candidate_telomere.fasta > ont_genome_tel.paf
## For hifi reads
minimap2 -t80 -x map-hifi ref_genome.fasta candidate_telomere.fasta > ont_genome_tel.paf
## Use paftools.js to process alignments in the PAF format. And use both Racon and Merfin algorithms to polish.
03.Identification_of_centromere_and_telomere_sequence
rf genome1.fa 2778010502000 -f -d -m -l 15
cd-hit -i genome1_filt.fa -o genome1_filt1.out -c 0.8 -T4 -n 4 -d 30
Use mafft to align the sequence and then ‘FastTree -nt -gamma -boot 500’to make the phylogenetic tree.
04.Genome_evaluation
1. BUSCO
busco -i genome.fasta -l embryophyta_odb10 -o busco -m genome --cpu 2
busco -i longest_isoform.fasta -l embryophyta_odb10 -o busco -m prot --cpu 2
2. N50 length
java -jar gnx.jar genome.fasta >genome.fasta.N50
3. qualimap2
bwa index genome.fasta
bwa mem -t 128 -M genome.fasta ${Sample}_1_clean.fq.gz ${Sample}_2_clean.fq.gz >${Sample}.sam
samtools view -@ 12 -bS ${Sample}.sam -o ${Sample}.bam
samtools sort -@ 12 ${Sample}.bam -o ${Sample}.sort.bam
qualimap bamqc -bam ${Sample}.sort.bam -nt 12 --java-mem-size=80G -c -outdir ./${Sample}
minimap2 -t 12 -ax asm20 genome.fasta ${Sample}.ccs.fasta > ${Sample}.minimap.sam
samtools view -@ 12 -bS ${Sample}.minimap.sam -o ${Sample}.minimap.bam
samtools sort -@ 12 ${Sample}.minimap.bam -o ${Sample}.minimap.sort.bam
qualimap bamqc -bam ${Sample}.minimap.sort.bam -nt 128 --java-mem-size=80G -c -outdir ./${Sample}
4. QV
meryl k=21 count output ${Sample}_1.meryl ${Sample}_1_clean.fq.gz threads=12
meryl k=21 count output ${Sample}_2.meryl ${Sample}_2_clean.fq.gz threads=12
meryl union-sum output ${Sample}.meryl ${Sample}_*.meryl threads=12
merqury.sh ${Sample}.meryl genome.fasta QV_${Sample}
5. LAI
gt suffixerator -db genome.fasta -indexname ${Sample} -tis -suf -lcp -des -ssp -sds -dna
gt ltrharvest -index ${Sample} -minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 1 -similar 85 -vic 10 -seed 20 -seqids yes > ${Sample}.harvest.scn
LTR_FINDER_parallel -seq genome.fasta -threads 30 -harvest_out -size 1000000
cat ${Sample}.harvest.scn genome.fa.finder.combine.scn > genome.fa.rawLTR.scn
LTR_retriever -genome genome.fasta -inharvest genome.fa.rawLTR.scn -threads 40
LAI -genome genome.fasta -intact genome.fa.pass.list -all genome.fa.out -t 20
05.Genome_annotation
1. TE annotation
perl EDTA.pl --genome genome.fasta --species others --overwrite 1 --sensitive 1 --anno 1 --evaluate 1 --threads 64
RepeatMasker -e rmblast -s -gff -nolow -no_is -norna -lib genome.fasta.mod.EDTA.TElib.fa -pa 64 genome.fasta
2. Protein-coding gene structure prediction
glimmerhmm genome.fasta -d trained_dir/watermelon -g -f -n 1 -o ${Sample}_glimmerhmm
augustus --strand=both --genemodel=complete --singlestrand=false --noInFrameStop=true --protein=on --introns=on --start=on --stop=on --cds=on --codingseq=on --alternatives-from-sampling=true --gff3=on --outfile=${Sample}.gff --uniqueGeneId=true --species=Watermelon genome.fasta
braker.pl --cores 48 --species=Watermelon --genome=genome.fasta --softmasking --workingdir=${Sample} --bam=${bam1},${bam2}... --gff3 --useexisting
exonerate -q ${file} -t genome --model protein2genome --querytype protein --targettype dna --showvulgar FALSE --softmaskquery yes --softmasktarget TRUE --minintron 30 --maxintron 10000 --showalignment FALSE --showtargetgff TRUE --geneseed 250 --score 800 --bestn 1 --verbose 0 --saturatethreshold 100 --dnahspthreshold 60 --dnawordlen 14 >${Sample}.gff
stringtie -p 30 -o ${Sample}.merge.gtf ${Sample}.merge.bam
gtf_genome_to_cdna_fasta.pl ${Sample}.merged.gtf genome.fasta > transcripts.fasta
gtf_to_alignment_gff3.pl ${Sample}.merged.gtf > transcripts.gff3
TransDecoder.LongOrfs -t transcripts.fasta
TransDecoder.Predict -t transcripts.fasta
cdna_alignment_orf_to_genome_orf.pl transcripts.fasta.transdecoder.gff3 transcripts.gff3 transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3
augustus_GFF3_to_EVM_GFF3.pl ${Sample}.augustus.gff > ${Sample}.augustus.gff3
glimmerHMM_to_GFF3.pl ${Sample}.glimmerhmm.gff >${Sample}.glimmerhmm.gff3
Exonerate_to_evm_gff3.pl ${Sample}.exonerate.gff >${Sample}.exonerate.gff3
partition_EVM_inputs.pl --genome genome.fasta --gene_predictions ${Sample}.ABINITIO_PREDICTION.gff3 \
--protein_alignments ${Sample}.exonerate.gff3 --transcript_alignments transcripts.fasta.transdecoder.genome.gff3 --segmentSize 100000 --overlapSize 10000 \
--partition_listing partitions_list.out
write_EVM_commands.pl --genome genome.fasta --weights `pwd`/weights.txt --gene_predictions ${Sample}.ABINITIO_PREDICTION.gff3 --protein_alignments ${Sample}.exonerate.gff3
--transcript_alignments transcripts.fasta.transdecoder.genome.gff3 --output_file_name evm.out --partitions partitions_list.out > commands.list
execute_EVM_commands.pl commands.list |tee run.log &
recombine_EVM_partial_outputs.pl --partitions partitions_list.out --output_file_name evm.out
convert_EVM_outputs_to_GFF3.pl --partitions partitions_list.out --output evm.out --genome genome.fasta
find . -regex ".*evm.out.gff3" -exec cat {} \; > EVM.all.gff3
3. Non-coding RNA gene annotation
tRNAscan-SE -o ${Sample}.tRNA.out -f ${Sample}.tRNA.ss -m ${Sample}.tRNA.stats -a ${Sample}.tRNA.fa genome.fasta --thread 20
rnammer -S euk -multi -f rRNA.fasta -h rRNA.hmmreport -xml rRNA.xml -gff rRNA.gff2 genome.fasta
cmscan -Z 100 --cut_ga --rfam --nohmmonly --tblout ${Sample}.tblout --fmt 2 --cpu 20 --clanin Rfam.clanin Rfam.cm genome.fasta > genome.cmscan
perl infernal-tblout2gff.pl --cmscan --fmt2 genome.tblout >genome.infernal.ncRNA.gff3
4. Gene function annotation
makeblastdb -in all.plant.protein.faa -dbtype prot -parse_seqids -title RefSeq_plant -out plant
makeblastdb -in uniprot_sprot.fasta -dbtype prot -title swiss_prot -parse_seqids out swiss_prot
blastp -query ${Sample}.protein.fa -out ${Sample}.swiss_prot.xml -db swiss_prot -evalue 1e-5 -outfmt 5 -num_threads 20
blastp -query ${Sample}.protein.fa -out ${Sample}.RefSeq_plant.xml -db plant -evalue 1e-5 -outfmt 5 -num_threads 20
interproscan.sh -f TSV,XML,GFF3 -goterms -cpu 1 -i ${Sample} -iprlookup -pa -T ${Sample}_temp -b ${Sample}
往期部分文章
「1. 最全WGCNA教程（替换数据即可出全部结果与图形）」
「2. 精美图形绘制教程」
「3. 转录组分析教程」
「4. 转录组下游分析」
「小杜的生信筆記」 ，主要发表或收录生物信息学教程，以及基于R分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料!
基因组
苏公网安备32011502012024号