1 材料与方法
1.1 样品采集与保存
注意要尽可能多的涵盖所有现生类群——所关注的阶元。
（eg：长蝽总科研究时的科级阶元）每个科中包含多个代表物种，亚科阶元中尽量分布均匀，含有属数量较多的亚科可以视情况多取。
重点关注系统发育位置存在争议或位置不明确的类群。
对样品进行定种，不能确定种的尽量确定属。（尽可能保存同一采集地点、同一形态形的同种样品作为二级凭证进行保存，二级凭证为干制标本或保存在>95%浓度的乙醇中。）
将标本保存在RNAlater中（RNAlater：样品比例 = 5：1）
物种列表S1中为什么有许多个体数目>1？
1.2 对已发布的转录组数据与基因组数据进行利用
2 核苷酸序列数据（转录组）
2.1 提取RNA
根据RNA提取试剂盒的说明，使用Trizol试剂提取各物种的总RNA
使用生物分析仪与分光光度计对提取出的RNA进行质检，如浓度、片段大小、RNA 完整性编号（RIN）、28S/18S 和 OD（260/280；260/230）等。
如果质量不合格是重新送样还是放弃使用该物种？
合格的 RNA 样品（RIN > 6 且总量 > 5 µg）将被用于文库构建
RIN（RNA完整值）是安捷伦公司开发的以数字化形式表示 RNA 完整性的参数，基于整体电泳图，将真核生物总RNA质量分类，并设定RIN数值范围为1-10，其中1代表降解最严重的情况，10代表最完整。
2.2 cDNA文库构建
根据试剂盒上的指示进行操作构建文库：
①在72℃下使用RNA裂解试剂将长链RNA剪切成短片段
②在裂解的RNA片段中加入逆转录酶使其合成第一链cDNA
③使用RNase H降解双链中的RNA
RNase H：催化 DNA - RNA 引物杂合体的RNA部分的降解的核糖核酸酶
④使用DNA聚合酶I合成第二链cDNA
⑤cDNA的甲基化*
⑥末端修复/加A尾
⑦添加接头
2.3 转录组测序
使用Agilent 2100 生物分析仪和 ABI StepOnePlus Real-Time PCR 仪对cDNA 片段的大小和浓度进行验证。cDNA 文库随后根据制造商的协议在 Illumina HiSeq2000 上进行测序，每个文库都进行了150bp双端测序，生成原始数据。
两端的序列是保护碱基（terminal sequence）、接头序列（adapter）、索引序列（index）、引物结合位点（Primer Binding Site）：其中 adapter是和flowcell上的接头互补配对结合的；index是一段特异序列，加入index是为了提高illumina测序仪的使用率，因为同一个泳道可能会测序多个样品，样品间的区分就是通过index区分。
Read：
对于每个通过质控参数的簇，一个序列被写入相应样本的R1 FASTQ文件，而对于双端测序运行，另外一个序列也被写入该样本的R2 FASTQ文件。 FASTQ文件中的每个条目包含4行：
- 序列标识符，其中包含有关测序运行和簇的信息。 该行的具体内容会因使用的BCL到FASTQ转换软件而不同。
- 序列（碱基信号； A，C，T，G和N）。
- 分隔符，只是一个加号（+）。
- 读取碱基的质量值。 这些是Phred +33编码的，使用ASCII字符表示数字质量值。
质量值Q是测序错误率的对数*-10。假如错误率是0.01，则Q值为20。可见，错误率越低，其Q值越高。即Q值越高越可靠。
0-31位都是控制字符，没法打印和保存，能打印的从要从32位的Space开始，所以就可以给实际的Q值加上一个固定值，这样就可以打印出来而保存在fastq文件中了。
所以下面就是固定值加多少的问题。phred33，就是加了33，也就是10变成43，查表是+。例子中很多F,Q值是70-33=37。（旧版的illumina加64位可使用脚本辅助人工来判断判断是Phred33还是Phred64，详见fastq格式文件及phred33的判断 - 简书 (jianshu.com)）
还是看不懂这个质量值是怎么搞得，迷迷糊糊
2.4从头组装转录组数据（核苷酸序列）
使用 SOAPdenovo-Trans-31kmer 1.01 组装转录本。
对污染数据进行预处理，去除质量较差的读数。
预处理
①去除有接头污染的reads（最小比对长度：15 bp；最多 3 个错配）
②N总数>10的reads
③大于 50 个碱基对的低质量reads（Phred 质量得分 = 2，ASCII 66 "B"，Illumina 1.5+ Phred+64）
其余数据用于从头开始的转录组组装（SOAPdenovo-Trans-31kmer，四步）
（1）作图前处理：
（进行kmer分析，用以了解：①基因组杂合度②重复序列比例③预估基因组大小④测序深度）
所有reads被分为31-mer以构建de Bruijn图；排除含Ns的K-mers；如果某个 K-mer有一个以上的出度，则剔除丰度小于 10%的出度。
（2）形成contigs（重叠群）
i 的入度（in-degree）是指向 i 的边的数量，出度（out-degree）是远离 i 的边的数量。
剔除丰度小于总出度的 5%或小于总入度的 2%的边。随后，边缘被报告为等位基因。
（3）Mapping
利用成对末端序列，将所有短读序列锚定到重叠群上，重建所有重叠群之间的关系。
（4）Scaffolding
根据它们之间的关系，使用所有大于 100 bp 的contigs构建scaffold。使用 Gapcloser（SOAPdenovo 2 软件包的一部分）填补scaffold中的所有gaps。
2.5 组装后转录组序列的质量评估
利用VecScreen和UniVec筛选载体和连接子/接头污染。
UniVec是一个数据库，可用于快速识别核酸序列中可能来自载体来源(载体污染)的片段，还包含用于克隆cDNA或基因组DNA过程中常用的adpter、linkers和引物的序列。
VecScreen是一个系统，它可以快速找到核酸序列的片段，这些片段可能来自于载体。
VecScreen在一个查询序列中搜索匹配UniVec中任何序列的段。
转录组 DNA 序列与已公布的载体 DNA 序列匹配并不一定意味着已确定的转录组 DNA 序列是污染物。Illumina 测序不需要克隆载体。因此，阳性结果可能是假阳性，也可能是外来遗传物质融入基因组的结果，如病毒基因组或细菌基因的一部分。虽然后者可以为进化基因组分析提供有价值的信息，但它会干扰转录本的同源分配。
因此，与载体或连接体 DNA 序列的末端强匹配或中度匹配被移除。在内部匹配的情况下，移除潜在的 DNA 序列污染，并将转录本分成两个独立的 DNA 序列。
所有至少 200 bp 的序列都被提交到 NCBI中的TSA数据库，以期进一步筛选外来污染序列。
（如何使用TSA数据库进行筛选？）
将含有污染物的序列剔除后剩余序列用于后续分析。
2.6 直系同源基因比对
使用基于profile hidden Markov models（pHMMs）算法的软件HaMStR将转录本以单拷贝形式映射到已知参考物种基因中
其中设置了78个参考类群，也可向其中添加新的。基因/基因列表。
使用 HMMER 软件利用 pHMM 在翻译水平上搜索符合 pHMM 的转录本，即蛋白序列。
随后用 BLAST 软件对照参照物种的全套基因搜索候选转录本，以剔除错误的转录本。
在OrthoDB 5 数据库（http://cegg.unige.ch/orthodb5）中检索需要比对类群的单拷贝直系同源基因以及相应氨基酸序列以及可用的基因描述和 InterPro注释信息。
然后使用 MAFFT的 L-INS-i 对齐算法对每个 OG 的氨基酸序列进行对齐。
然后，以得到的多氨基酸序列比对结果为基础，使用 HMMER 3.0 软件包中的 hmmbuild 程序构建 pHMM。
为了对候选转录本进行反向 BLAST 搜索，我们使用 12 个参考物种的官方基因组生成 BLAST 数据库，将这些基因组与 OrthoDB 5 分析的基因组相对应。
HaMStRad进行转录本的直同源评估（省略）
使用宽松的核苷酸比对检索标准……
2.7 同源基因的功能注释
通过基因本体论（Gene Ontology，GO）术语对 1,478 个同源蛋白编码基因的可能功能进行描述。
GO总共有三个ontology（本体），分别描述基因的分子功能（molecular function）、细胞组分（cellular component）、参与的生物过程（biological process）。GO的基本单位是term（词条、节点），每个term都对应一个属性。
首先使用 Blast2GO 2.5.1为 12 个参考物种官方基因组的所有基因分配GO 术语。
12个参考物种：
Ixodes scapularis(Arachnida: Ixodida), Daphnia pulex (Branchiopoda: Diplostraca）Zootermopsis nevadensis (Insecta: Isoptera), Pediculus humanus (Insecta: Phthiraptera）Acyrthosiphon pisum (Insecta: Sternorrhyncha）Nasonia vitripennis, Apis mellifera, Acromyrmex echinatior (Insecta: Hymenoptera), Tribolium castaneum (Insecta: Coleoptera), Bombyx mori (Insecta: Lepidoptera), Anopheles gambiae and Drosophila melanogaster (Insecta: Diptera).
如果是无参基因组的物种，就需要根据unigene的Nr注释，用Blast2go软件得到unigene的GO注释信息。
利用 Blast2GO，我们对官方基因组的所有氨基酸序列与 NCBI 的非冗余氨基酸序列数据库（nr）（设置：E 值：1E-03；HSP 长度截止值：33）进行了 BLAST相似性搜索（使用 BLASTp）。（目的：在 nr 数据库中为官方基因组的每个序列找出最相似的序列。所有确定的序列都与附有证据代码的 GO 术语相关联。因此，只有 BLAST 找到的序列才能用 GO 术语进行注释。）
从整体上来看，GO注释系统是一个有向无环图(Directed Acyclic Graphs)，GO各term之间的关系是单向的，GO term之间的分类关系有三种：is a、part of 和 regulates。
其中，Pvalue 这一行，如果大于0.05，即会显示NA，即图中只显示显著的P value。
形状的含义：程序默认把显著性最高的前10个GO term设置为方形，其他的GO term为圆形。
颜色的含义：颜色越深，代表该GO term越显著。颜色由浅到深分别为：无色——浅黄——深黄——红色。
使用 R 统计软件 2.15.1 版的 TopGO 软件包中的消除算法（elim algorithm）和费舍尔精确检验（Fisher's exact test）比较 12 个参考物种中每个物种的直系基因集（1,478 个直系基因子集）与相应官方基因集的序列，确定了分配给 GO 领域 "生物过程 "的高比例 GO 术语。
2.7 数据沉积
所有分类群都作为样本被保存在NCBI 生物样本数据库中；对于在 NCBI 分类学中未知的物种也进行了注册，并获得了一个新的分类标识（Tax-ID）；原始成对末端读数和组装后的转录组已提交至序列读数档案（SRA）和NCBI的转录组枪组装（TSA）数据库。
所有数据（包括实验描述）均可通过 NCBI 的Umbrella BioProject ID PRJNA183205获取
4617
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
