转录组原理

一个项目测多少数据量合适
一个样品(sample)测6-10G数据
6G → 碱基的数目，双末端测序对应20M数据（20M为reads条数）20M x 300 (PE每条数据300个碱基，左右各150)
国内说测数据量，以G为单位；发文章标准一般说测多少M数据，不说G(涉及到测序原理，抽样以DNA片段为样本抽样)
当数据量大于20M后，转录本数量基本饱和，再增加数据量转录本增加数量很少（意义不大，性价比不高）
若研究mRNA则没有必要很高的数据量，10G就太高了
饱和度曲线

饱和度曲线

转录组分析步骤

0 准备基因组序列和测序数据

-1 比对把各个样本的fastq格式的reads比对到基因组序列上，得到一个bam格式的文件（sample.bam）

软件：Hisat2(动植物推荐，精确性稍差，速度快很多)、STAR(最好的软件，运行内存要求高，精确度高，医学重用)
得到一张表格，比对率表格，大部分动植物比对率至少应该在70%以上
比对到转录组，用bowtie2, 不推荐使用，已经不是标准做法
conda install hisat2
比对步骤：hisat2-build 构建参考基因组（构建到参考基因组文件夹下）
→ hisat2比对(得到一个 .sam 格式的文件)
→ sam2bam(文件压缩、排序过程) 命令 samtools sort
→ samtools index (bamindex构建index索引) 命令 samtools index
注： hisat2-build 过程耗时长(大概要十几个小时)，消耗内存大(像人类基因组用STAR就需要消耗至少16G内存)
有些物种不需要构建参考基因组，在hisat2官方网站右边已列出已经建好的参考基因组，可以直接下载。
by the way, 官网上人类参考基因组有：genome / genome_snp(结合了SNP构建的参考基因组) / genome_tran(结合了一些已知的基因结构构建的参考基因组) / genome_snp_tran(前两个都有) 建议选择最后一个
hisat2文件挂起nohup有日志文件输出，后面加1>xxx.log 2>xxx.err表示将运行过程的标准输出输出到xxx.log文件，将错误输出到xxx.err文件中；后面加1>xxx.log 2>&1表示将错误输出也输入到xxx.log文件
bam格式文件为二进制文件(binary file), cat、less等都打不开，可以通过 samtools view xxx.bam | less -S 查看 (sam文件中为乱序，bam文件按染色体序列排序)
bamindex那一步需要装 IGV —— 对比对结果可视化的软件
软件分类：①基于比对的软件：feautre count; ②基于Kmer频率的软件：Salmon / Killing stone 都可以用，后者较新，建议用前者，因为文章中都需要比对率，后者不给出比对率。

-2 定量数每个基因上落了几条reads，需要将基因结构画在染色体上 → 基因表达量（表格），又称原始的reads count 矩阵
以上为转录组标准分析
（非模式物种的可变剪接、融合基因分析不值当）

g/S	Sample1	Sample2	Sample3
gene1	38	55	76
gene2	127	41	86
gene3	46	29	34
...	...	...	...

在基因层面定量，还是在isoform层面（剪切形式）做定量
软件：htseq → htseq-count / Subread → featureCounts（推荐，新出）conda 安装
conda install bioconductor-rsubread
需要输入的文件有sample.bam(描述每条序列比对到的基因组的位置)、genes.gtf(基因组注释文件，描述每个基因的起始、终止位置)
输出文件 → 没有固定格式，表达量文件 feature.count
用subread软件表里的featurecounts程序

stringtie 用于做转录组拼接，需要证明自己的拼接靠谱，一般不用，除非证明拼接靠谱，已经发表并被广泛认可

-3 原始reads count矩阵标准化（细节见转录组原理篇课程）
RPKM(PE)/FPKM(SE) （该方法是错的，但得出的结论基本上是对的）
TPM (对的)
(重复序列处理方式：可以匹配到多个位置。uniq, EM算法，自动用贝叶斯算一个概率)

最后编辑于：2020.07.14 08:07:45