文献笔记四十三：不同形态的南瓜重测序探索与形态和有价值的农艺性状有关的基因组变异

文章题目

Whole-genome resequencing of Cucurbita pepo morphotypes to discover genomic variants associated with morphology and horticulturally valuable traits

期刊、发表时间、完成单位

Horticulture Research 2019年8月 University of Thessaloniki Greece 当前阶段重点关注基因组重测序数据分析方法以及对结果的解释

试验材料

图片来自原论文

8个品种 为了整合品种内变异（To integrate putative intra-cultivar variation），每个品种采集5株植株的幼叶然后混合。使用SSR引物评估品种内的变异。

DNA测序

• 测序平台：BGISEQ-500
• 双端150bp

原始数据处理、与参考基因组比对、snp检测

• 去除接头（The adapter sequence in the raw data was removed.）
• 低质量碱基或者含有太多的未知碱基的reads去掉
• 与参考基因组比对使用BWA软件
• SNP检测使用SOAPsnp SAMtools 和 GATK三个软件
• AnnoDB软件注释变异结果
• 对结果进行过滤 深度小于10，质量值小于20，缺失数据小于30%，非二等位基因

利用得到的snp结果计算一些参数，包括：

• 遗传多样性 Genetic diversity TASSEL
• 连锁不平衡 linkage disequilibrium Plink v1.9
• snp数量 https://github.com/JoseBlanca/variation/
• 杂合位点和纯合位点的比例 https://github.com/JoseBlanca/variation/

以上每个品种单独计算，然后分亚种计算（number of SNPs were calculated for the eight accessions, as well as for each subspecies separately.） (关于连锁不平衡的一系列计算还有点看不明白)

• Genetic pairwise differentiation
• 无根网络 SplitsTree
• Weir and Cockerham's Fst 表示亚种间的遗传距离 vcftools

候选基因的遗传变异分析

• 使用snpEff软件注释结果
• 根据之前的研究结果挑选感兴趣的数量性状位点
• 37个与开花，果实果肉颜色，叶片形态有关的候选基因
• 18个与果实形状有关的候选基因
• 近来的研究表明携带TRM(Tonneau1 recruiting motif)的蛋白与OVATE相互作用在子房发育过程中调节细胞分裂模式从而改变最终的果实形状。
• blast方法鉴定带有TRM结构域的蛋白，evalue设置10-50
• 序列比对IQ-tree构建进化树 两种方法评估树的支持率
• 鉴定经受不同选择压力的候选基因 计算Fst和Tajima's D 使用vcftools

试着重复其中的一些分析结果

下载参考基因组和原始测序数据

在论文中可以找到对应的下载链接

• 数据量相对来说还是比较大的，我只选取原始测序数据的前4000000行；参考基因组只选择第一条染色体 从基因组文件中挑出一条染色体使用的是samtools命令

samtools faidx Cpepp_genome_v4.1.fa
samtools faidx Cpepp_genome_v4.1.fa Cp4.1LG01 > OneChrom_Cpepp_genome.fa

samtools faidx OneChrom_Cpepp_genome.fa
bwa index OneChrom_Cpepp_genome.fa

• 比对生成 bam文件

for i in Acorn Cocozelle Crookneck Marrow Pumpkin Scallop Yellow_Zuchini Zuchini
        do
                bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:'${i} Reference/OneChrom_Cpepp_genome.fa ${i}/${i}_1.fastq ${i}/${i}_2.fastq | samtools v
        done

这里有一个疑问：bwa mem 比对需要制定测序平台，原文测序平台BGISEQ-500，不知道这里设置为illumina是否可以？

• 为参考序列生成一个.dict文件

gatk CreateSequenceDictionary -R  OneChrom_Cpepp_genome.fa -O OneChrom_Cpepp_genome.dict

• 生成gvcf文件

for i in Acorn Cocozelle Crooknect Marrow Pumpkin Scallop Yellow_Zuchini Zuchini
        do
                samtools sort -@ 4 -m 4G -O bam -o Output/${i}.sorted.bam Output/${i}.bam
                echo 'sorted done'
                gatk MarkDuplicates -I Output/${i}.sorted.bam -O Output/${i}.sorted.markdup.bam -M Output/${i}.sorted.markdup_metrics.txt
                echo 'MarkDup done'
                samtools index Output/${i}.sorted.markdup.bam
                echo 'index done'
                gatk HaplotypeCaller -R Reference/OneChrom_Cpepp_genome.fa --emit-ref-confidence GVCF -I Output/${i}.sorted.markdup.bam -O GVCFoutput/$
                echo 'Haplo done'
        done

• 合并g.vcf文件 这里参考GATK 4.0 全外显子call variant

gatk CombineGVCFs -R ../Reference/OneChrom_Cpepp_genome.fa $(for i in $(ls *.vcf);do echo "--variant $i"; done) -O Combined.g.vc

• 检测变异

gatk GenotypeGVCFs -R ../Reference/OneChrom_Cpepp_genome.fa -V Combined.g.vcf -O Nangua8cultivars.vcf

• 查看生成的vcf文件中有多少个变异

bcftools view -H Nangua8cultivars.vcf | wc -l

结果为 224755

以上步骤就得到了计算一些群体基因组学相关指标所用到的vcf文件

这里为了减小运算压力，只选取了原始数据的前400000行和参考基因组的第一条染色体

没有对原始数据进行质控过滤

今天就先到这里啦！