测序数据质量控制
测序数据分析前需要经过数据预处理，并检查数据GC含量、序列重复、是否存在接头等。
- 质量评估：
使用 FastQC 检测原始数据质量
fastqc –o fastqc_results –f fastq test_1.fastq test_2.fastq b_1.fastq b_2.fastq
GC含量分布图.问题区域需要切除，例如切除前 30bp
错误率分布图.越接近尾部错误率越高，通常要求 Quality Score > 30
- 质量控制
使用 Trimmomatic 去除低质量reads。
Trimmomatic 详细说明参考：https://www.jianshu.com/p/a8935adebaae
FastQC和Trimmomatic的安装及使用参考：https://www.jianshu.com/p/bc3ad9379e3e?utm_campaign=hugo&utm_content=note&utm_medium=seo_notes&utm_source=recommendation
用法：
java -trimmomatic.jar PE -threads 2 -phred33 \test_1.fq.gz test_2.fq.gz \test_1.trimed.fq.gz test_1.un.fq.gz test_2.trimed.fq.gz test_2.un.fq.gz \ILLUMINACLIP:/path/to/Trimmomatic/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:76
在质控后，再质检一次，对比看看有什么不同。
reads比对
将 reads 匹配到参考基因组或转录组的相应位置上
• 非剪接比对：转录组
Bowtie、BWA
• 剪接比对：参考基因组
STAR、HISAT、Topha
对鉴定SNP做了优化：GSNAP、MapSplice等
HISAT2比对流程
① 建立基因组索引
extract_splice_sites.py tair10.gtf > genome.ss # 把剪切位点提取出来extract_exons.py genome.gtf > genome.exon # 把exon提取出来hisat2-build --ss genome.ss --exon genome.exon genome.fasta genome # 最后的genome是输出文件的前缀
②利用注释文件比对
hisat2 -p 4 --known-splicesite-infile genome.ss --dta -x tair10 -1 test_1.trimed.fq.gz -2 test_2.trimed.fq.gz -S test.sam ## -p 线程数 ## --known-splicesite-infile 输入剪切位点文件## --dat 转录本拼接## -x index 库文件前缀CDS 和 exon 前## -1 -2 双端测序 fastq的名字， 如是单端测试 –U ## -S 输出文件，是比对的 SAM 文件
没有注释文件的比对方法：
hisat2 -p 18 --dta -x ~/genome/rice -1 /path/to/Rice_1.fq.gz -2 /path/to/Rice_.fq.gz -S rice.sam
③ SAM 文件处理
使用 samtools 对 SAM 文件排序并转化为 BAM 文件。samtools是一个用于操作sam和bam文件的工具合集，包含有许多命令。
samtools view -bS SRAxxx.sam > SRAxxx.bam # 查看bam文件内容samtools sort -@ 2 -o SRAxxx.sort.bam SRAxxx.bam # 按比对位置排序+格式转换samtools index rice.bam # 建立bam文件索引samtools merge -@ 4 -h SRR1582649.bam merged.bam SRRxxx1.bam SRRxxx2.bam SRRxxx3.bam # 把生成的bam文件合并为一个文件## -@ 额外线程数## -m 每个线程最大占用内存，单位 K/M/G，根据实际情况调整。## -o 输出文件## -h 因为每个文件的sam文件表头都一样，所以用-h指定某一个文件的表头作为总文件的表头
④比对结果可视化
比对结果使用 IGV 、Genome Maps 和Sacant 等可视化查看。
例如：IGV 通过读入基因组和注释信息以及BAM 文件展示比对结果。
需要额外添加 BMA 的索引:samtools index、test_sorted.bam、test_sorted.bai
⑤比对结果评估
比对结果评估工具：RSeQC、Qualimap
- Reads 匹配百分比评估预测精度和DNA污染程度或参考基因组的选择是否适合；
- Reads 随机性分布 评估reads打断的随机程度；
- 匹配Reads的GC含量，与PCR偏差有关。
基于NGS的转录本定量---StringTie
- reads 计算策略
① 只选唯一匹配 reads：用于估计基因水平的 reads 匹配数，常用工具如HTSeq-count、featureCounts；
② 保留多重匹配的 reads：利用统计算法将多重比对reads定位到对于的转录本异构体上，如 Cufflinks、StringTie、RSEM等
用法：
stringtie merge.bam -e -G Oryza_sativa.genomic.gff -o Oryza_sativa.new.gtf
RNA定量标准化
Reads 数受到基因长度、测序深度和测序误差等影响，需要归一化处理。
常用归一化策略：
RPKM（reads per kilobase of exon model per million reads）每 100 万 reads 比对到每 1kb 碱基外显子的 reads 数目
FPKM （fragments per kilobase of exon model per million reads） 每 100 万 reads 比对到每 1kb 碱基外显子的片段数目
TPM（transcripts per million）每 100 万 reads 比对到该转录本的片段数目
单端测序：RPKM = FPKM
双端测序：FPKM 更优
（FPKM 和 TMP 可以换算）
转录组核心分析策略选择
有参考基因组
可基于参考基因组序列和注释文件（gff或gtf文件）获得注释基因的基因表达、可变剪接状况、识别新转录本；
• 分析策略：
① 利用 TopHat、HISAT 及 STAR 等（剪接匹配工具）将 reads 比对到参考基因组上；
② 使用 Cufflinks 和 Stringtie 等软件进行转录本拼接，获取转录本表达量；
③ 拼接结果还可以与注释文件相比获取新转录本，并对新转录本进行功能注释。
无参考基因组
以前的策略：是对二代测序 reads 使用 Tiinity、SOAPdenovo-Trans等拼接成转录本，再获取表达信息。
• 现在的策略：
① 使用三代测序（Pacbio）进行全长转录本测定，获取参考转录组；
② 对获得的全长转录组进行功能注释；
③ 进行二代测序，并利用 Bowtie 等非剪接比对工具直接匹配到到参考转录组上；
④ 使用 RSEM、HTSeq-count 等工具进行转录本定量。
以 stringtie 进行有参转录本拼接
① 转录本拼接
stringtie -p 4 -G reff.gtf -o test.gtf -l test test_sorted.bam## -p 线程数## -G 参考基因组注释文件## -o 输出文件## -l 转录本命名前缀
② 转录本整合
stringtie --merge -p 4 -G ref.gtf -o test_merge.gtf gtflist.txt## 先将不同样本在上一步产生的gtf文件路径写进 gtflist.txt，再合并 gtf，用于后续分析鉴定
③ 转录本注释文件比较
## 先将不同样本在上一步产生的gtf文件路径写进 gtflist.txt，再合并 gtf，用于后续分析鉴定gffcompare –r ref.gtf -G -omerged test_merge.gtf
FPKM的计算
使用 stringtie 计算基因和转录本的 FPKM
stringtie -e -p 4 -G test_merge.gtf -A test_merge_gene.gtf -o test_merge_transcript.gtf test test_sorted.bam## -e 只列出已知转录本丰度，不进行额外拼接预测## -G 如果只关注参考基因组注释基因就用参考基因组gtf，如果关注新转录本就用前一步meregd的gtf## -A 输出有 gene 表达量的 gtf## -o 输出有转录本表达量的gtf
Reads 计数
## 使用 HTSeq-count 从比对结果中提取所有基因匹配的reads数据，用于后续差异表达分析htseq-count -q -f bam -s no -i gene_name /path/to/alignment/test_sorted.bam /path/to/genome/gff/tai /path/to/ref.gtf > test.count## -q 不显示进程报告## -f 比对文件格式 sam / bam## -s 是否考虑链特异性## -i 提取属性名
Q&A
- 定量过程中，比对到基因组多个位置的reads（multi mapped）是如何处理的？
A: 这些 reads 比对到基因组多个位置，无法确定是由哪个基因转录而来，所
以这部分 reads 在定量过程中直接被过滤掉。
- 有重叠区域的两个基因，重叠区域的reads在定量时如何分配？
A: 基因重叠区域的 reads 也无法确定由哪个基因转录而来，这部分的 reads
也是被过滤处理。
- 基因表达的阈值是多少？为什么设置为这个阈值？
A: 一般认为 FPKM 大于 1 时基因表达，这个阈值是主流杂志推荐的。
