首页 > 编程知识 正文

转录组上游分析和下游分析,转录组测序数据分析

时间:2023-05-05 00:27:52 阅读:179034 作者:4370

转录组分析目前常用于功能基因组领域。 相关软件组合种类也非常丰富,本文采用hisat2 samtools stringtie策略从转录组数据中挖掘差异表达基因。 现小编整理了该组合的运行流程,供日后查阅; 如果能在平台上分享,帮助更多的初学者,编辑将不胜荣幸。 如果有纰漏,希望各路大佬指正。

首先,让我们从整体上看一下软件执行的功能。

hisat2 :建立参考基因组索引,比对读我

samtools:sam2bam的转换

stringtie :估计转录本表达量

使用的数据结构如下,这里使用了酵母的转录组数据和参考基因组:

双端测序数据已经进行了快速QC滤波,本文没有涉及具体的滤波过程。 示例数据仅供参考。执行时要关注文件格式转化,以及各种格式下包含的生物信息。

@ bio cloud :~~/1223/ngs 2022 $ tree .gene _ data.CSV//差异表达基因样本相当于最终输出的标准答案。 (-- genome (-- yeast.fa /参考基因组文档(--yeast.gff//参考基因组注释文档(-yeast _ transcript ome.fa /) 通常,下机时公司会同时推出rawdata和cleandata,但直接使用后者就可以了。 (S1 _ y _1. FQ.gz )-- S1 _ y _2. FQ.gz-- S2 _ y _1. FQ.gz-- S2 ((((() ) () () ) ) ) ) ) )。 中可用的软件安装软件包fast QC _ v0. 11.9.ziphisa T2-2.2.1-Linux _ x86 _ 64.ziphisa T2-2.2.1-OS x

$ CD./genome $ hisat2-build yeast.fa genome.fa//创建参考基因组索引,且hisa T2-build之间不要有空格。 //在新创建的文件夹中进行匹配操作$ CD . $ mkdir alignment $ CD alignment $ hisa T2-P6-x ./genome/genome.fa-1 ./reads/S1 _ S1 $ samtoolsview-bs S1.Sam-os1.bam//samtoolsview将Sam文件转换为bam文件。 bam是sam的二进制文件,二进制转换后的存储容量会更小。 $ samtoolssorts1. ba ms1.sorted//samtoolssort对S1.bam进行排序,并生成s1.sorted.bam。 后者会自动添加. bam后缀,不需要将其添加到命令行中。 $ samtools index s1.sorted.bam //在排序的bam文件中输入$ string tie./S1.sorted.bam-g ./genome/索引yeast.gff-e-p2-o.string tie根据gtf注释和分级的bam的比较结果推测转录本表达量//用同样的方法得到了其他处理下的双端序列的比较结果: S2 这里不一一重复。 //在新创建的文件夹中手动使用$ CD . $ mkdir differential _ expression $ CD differential _ expression $ vim sample _ list.txt # //sample_list.txt文件如此长,这里我们将两个gtf文件放在了differential_expression文件夹下。 //S1/mnt/1223/ngs 2021/differential _ expression/S1 _ out.GTF//S2/mnt/1223/ngs 2021/differential _ ed

f$ python prepDE.py3 -i sample_list.txt -g gene_count.csv -t transcript.csv //stringtie的脚本生成差异表达基因列表。// 注意stringtie差异基因汇总脚本有prepDE.py和prepDE.py3两个版本,前者适用于python2环境,后者python3环境。使用混了会报错。本文基于python3.9.5环境。

执行完上述步骤后得到的gene_count.csv即为最终结果,可以导入到R语言用edgeR / DESeq2包进行差异表达基因分析。有机会再整理后继教程。

sam文件基础

sam文件在序列分析中至关重要,无论是转录组分析还是基因组call SNP。认识sam文件所包含的信息有助于理解数据。我们依然以本文中生成sam文件的命令行举例。转录组分析核心比对步骤是使用hisat2将测序得到的reads比对到建立好索引的参考基因组:

$ hisat2 -p 6 -x ../genome/genome.fa -1 ../reads/s1_y_1.fq.gz -2 ../reads/s1_y_2.fq.gz -S ./s1.sam -x 参考基因组 .fa文件-1 双端测序第1段 .fq.gz格式-2 双端测序第2段 .fq.gz格式-S 输出文件地址+名字,输出结果为SAM格式
双端测序一般为read1.fq / read1.fq.gz / read1.fq.bz2格式,前后两端同时出现,并同时作为hisat2输入。将测序得到的fastq文件经过hisat2比对到参考基因组得到SAM文件。SAM文件就是序列比对文件,有上下两个部分,分别包括头部注释部分和比对结果部分: 头部注释部分 //头部注释部分以@开头:@HD VN:1.0 SO:unsorted //HD行:VN的版本以及比对排序类型。此处显示SO:unsorted表示没有排列顺序。@SQ SN:NC_001133.9 LN:230218//SQ行:参考序列目录。SN:参考序列名字。LN:参考序列长度@SQ SN:NC_001148.4 LN:948066@SQ SN:NC_001224.1 LN:85779@SQ SN:NC_001140.6 LN:562643@SQ SN:NC_001141.2 LN:439888@SQ SN:NC_001142.9 LN:745751@SQ SN:NC_001143.9 LN:666816@PG ID:hisat2 PN:hisat2 VN:2.0.5 CL:"/mnt/bai/public/单薄的小松鼠/hisat2-2.0.5/hisat2-align-s --wrapper basic-0 -p 8 -x ../genome/genome.fa --dta -S ./s1.sam -1 /tmp/1909596.inpipe1 -2 /tmp/1909596.inpipe2"//PG行:使用的比对程序名,此处为hisat2 比对结果部分

比对结果部分每行表示一个read和参考基因组的比对结果信息。前11列为主列,包含了大部分重要信息。

SRR5511068.3 99 NC_001147.6 1002516 60 75M = 1002624 184 GTACTTAACATTCTTCTAATCATGTTAAAAGGTAAAACCTGGCCCATTTTACGATCGATCTGTAAAATCTTATAC AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEEAEEEEEEEEEEEAEEEEEEEAEEEEEEEEEEA AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YS:i:0 YT:Z:CP NH:i:1SRR5511068.3 147 NC_001147.6 1002624 60 76M = 1002516 -184 GATAAGTAAGCAATGGTGGTAATTGCAATATTTTGCATATGTGCACGAGAAGAACTATTTTGAAGTAAGATCACTG /EEEE<EAEEEEEE6E/EEEEEEEAEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEAAAAA AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:76 YS:i:0 YT:Z:CP NH:i:1SRR5511068.9 83 NC_001146.8 299462 60 75M1S = 299435 -103 ATTGGAAAAGAAAGTCGCGGCAAAGAGAAATGCCAATAAGACCGGGAATCAAAATTCTAAAAAGAAGAGTCAGAAG <EAAA6<A<E/EEA/A</AE/EEEAA<EEAE6AEEAAEEAEEEAAAEE/EEEEEEEEEEEE//EEEEEE/EAAAAA AS:i:-1 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YS:i:0 YT:Z:CP NH:i:1

主体部分以tab分隔从左到右依次为:

QNAME:Query Name,测序reads名称。如果是双端测序则会出现二次,两端各比对上一次。FLAG:数值为2的次方数或者其加和,每个2的次方代表一种情况。详情可以参考CSDN大佬的博客:https://blog.csdn.net/genome_denovo/article/details/78712972
或者简书大佬:
https://www.jianshu.com/p/ab133ee9712cRENAME:Reference Name,双端测序R1比对上的参考序列名,没比上就是*。POS:position,双端测序R1起始比对上的位置序号。MAPQ:Mapping Quality,质量分数,越高代表越准确。CIGAR:比对结果,M代表完全匹配。RNEXT:双端测序R2端比对情况,比对不上用*号,比对到同一段用=号。MPOS:双端测序R2端比对位置。ISIZE:文库插入长度,R2端位置-R1端位置+CIGAR处的值。SEQ:序列信息。QUAL:reads质量,用ASCII码表示。 最后总结一下:

转录组数据分析步骤不仅仅是比对产生差异表达基因,后继还会涉及差异表达基因的显著性分析,相关性分析,GO和KEGG分析等等。本文只是总结了转录组分析的上游步骤。而就上游步骤而言,可用的比对软件种类也非常多,本文只是借助hisat2+samtools+stringtie的经典步骤来学习转录组上游分析。其他比对软件日后如有涉及再行整理。

版权声明:该文观点仅代表作者本人。处理文章:请发送邮件至 三1五14八八95#扣扣.com 举报,一经查实,本站将立刻删除。