一、 首先放一段illumina测序平台的官方原理链接：
https://www.bilibili.com/video/BV1Cx411p7dm
illumina测序平台采用的是二代测序原理，即"Next-generation"sequencing（NGS）其关键的特点是：高通量、边合成边测序
二、 illumina测序平台的主要步骤可大致分为4步：
1. 建库
A. 首先，进行3`末端突出的碱基A，为了后续可以在粘性末端添加接头（Adaptor）
B. 接着，利用突出的A-tail，添加环状结构以增强稳定性。
C. 然后，环状结构删除从而形成“Y”形接头。目的是为后续PCR中作为引物扩增继续添加文库index、与测序平台互补的寡核苷酸序列（作为后期测序引物）
D. 最后，利用不同PCR的扩增体系，将P7、P5以及Barcode加入到测序DNA上
Figure1
建库完成后的每条DNA的单链均一端连有测序引物Read1 Sequencing Primer（Rd1SP）和P5，另一端为Rd2 SP、Index（Barcode）和P7。
2. 桥式PCR扩增
“桥”式扩增（通过分别以P5和P7固定接头不断的扩增），将散布在表面的单核苷酸序列变成散布的DNA簇，以解决单分子产生的光信号很弱的缺点。
Figure2
3. 边合成边测序
在流通池加入可逆终止荧光dNTP，其3'-OH被阻隔（糖基3'连接有叠氮基团，在链延伸时起到了阻止添加下一个dNTP作用）。也就是说，达到一个循环只能扩增一个碱基，光信号捕捉仪只需要在一个信号结束的时候捕捉光信号，就可以确定该碱基。
Figure3
4. DNA簇颜色信号判断碱基
Figure4
三、 下机文件fastq的格式
一条序列有四行：
- 序列综合信息
- 序列碱基信息
- +号
- 序列对应的每个碱基的质量信息（根据荧光信号的强弱和弥散程度分析得到）
PS：phred33和phred64是针对read quality scores的转码形式命名的
Figure5
Figure6
Figure7
四、转录组数据的分析流程
以下步骤中所有软件的安装均基于linux系统，采用conda安装：
conda install fastqc
我更喜欢分环境安装不同的生信软件：
#将软件安装在特定环境下
conda install -n fastqc_env fastqc
#激活环境
conda activate fastqc_env
#检查是否安装成功
fastqc -v
#关闭环境
conda deactivate
主要用到的软件有：
- 数据质量分析软件 fastqc multiqc
- 数据质量控制软件 fastp(采用c++编写，功能强大，运行快)
- 序列比对软件 hisat2（以及samtools将sam文件转化为bam文件）
- 样本定量软件 stringtie
- 差异表达分析R语言包 deseq2 package或 ballgown package
以下代码以单个样本为例：
1. 测序结果Raw data的分析
#生成文件的md5值
md5sum *.gz > md5.txt
#比对md5值
md5sum -c md5.txt
#fastqc查看样本质量
fastqc -t 6 *gz -o result/
#multiqc整合fastqc分析结果
multiqc ./
# -t参数指定线程数 *gz为所有包含gz的文件 -o参数输出结果
#./指出于当前目录
2. Raw data的质量控制（转成Clean data）
#按照fastp默认参数进行质量控制，想增加参数可以查看官网maunal
fastp -i in.R1.fq.gz -I in.R2.fq.gz -o our.R1.fq.gz -O out.R2.fq.gz
# -i -I分别为双端测序的两个文件 -o为输出结果文件
3. 比对到物种参考基因组（获得bam文件）
#以羊的基因组为例，首先建立参考基因组索引
hisat2_extract_exons.py sheep.gtf > exons.txt
hisat2_extract_splice_sites.py sheep.gtf > ss.txt
#该命令运行所需内存很大超过200G
hisat2-build -p 6 --ss ss.txt --exon exons.txt sheep.fa sheep
#电脑配置不够可以采用精简命令（8G足矣）：
hisat2-build -p 6 sheep.fa sheep
#当使用精简命令后，在比对时要加入--known-splicesite-infile ss.txt参数
hisat2 --known-splicesite-infile ss.txt --dta --thread 6 -x sheep_ovi/sheep_ovi -1 1C.R1.fastq.gz -2 1C.R2.fastq.gz -S result/1C.sam
#将sam文件转化为二进制的bam文件，节约空间和处理时间
samtools view --threads 6 -m 2G -bS 1C.sam > 1C.bam
#将bam文件排序（必须进行这一步）：
samtools sort --threads 5 -m 2G 1C.bam -o 1C.sorted.bam
4. 样本定量分析（得到gtf文件）
第一种情况：
#该命令会影响后续差异分析结果，会出现很多STRING.xxx
stringtie -p 6 -G sheep.gtf -o 1C.gtf 1C.sorted.bam
stringtie --merge -p 6 -G sheep.gtf -o black_stringtie_merged.gtf mergelist.txt
stringtie -e -B -p 6 -G black_stringtie_merged.gtf -o sec_gtf/sample_1C/str_1C.gtf 1C.sorted.bam
第二种情况：
#不经过融合多个样本的结果这一步骤后，可以很好的避免差异分析中大量的STRING.xxx
stringtie -e -B -p 6 -G sheep.gtf -o third_gtf/th_1C.gtf 1C.sorted.bam
#后续用deseq2分析则需要该命令转化成deseq2可输入文件；用ballgown可忽略
prepDE.py -i deseq2_pre.txt -g black_deseq2_gene.txt -t black_deseq2_trans.txt
5. 分组差异表达分析
deseq2包的差异表达分析：
#读取文件
input_data <- read.table("black_deseq2_gene.txt", header = TRUE, row.names=1)
#取整
input_data <- round(input_data, digits = 0)
#准备工作
input_data <- as.matrix(input_data)
#指定分组名称
condition <- factor(c(rep("black",6),rep("kazakstan",6)))
#建立数据框
coldata <- data.frame(row.names = colnames(input_data), condition)
library(DESeq2)
#构建deseq输入矩阵
dds <- DESeqDataSetFromMatrix(countData = input_data, colData = coldata, design = ~condition)
#DEseq2进行差异分析
dds <- DESeq(dds)
#提取数据
res <- results(dds)
summary(res)
res <- res[order(res$padj), ]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds,normalized=TRUE)), by="row.names", sort=FALSE)
names(resdata)[1] <-"Gene"
#head(resdata)
#输出结果文件
write.table(resdata, file="diffexpr-results.txt", sep ="\t", quote = F, row.names = F)
ballgown包的差异表达分析：
#安装ballgown包
BiocManager::install("ballgown")
#载入包
library("ballgown")
#查看帮助
help("ballgown")
#引入文件
data_directory_pitu = file.path("C:/other_disk/E_disk_database/RStudiorun/rna_seq_bla_ovi/differ_gene/different_gene_ballgown/diffr_gene_ballg/extdata_pituitary")
# make the ballgown object：
bg_pitu = ballgown(dataDir=data_directory_pitu, samplePattern='sample', meas='all')
bg
#设置样本分组
sampleNames(bg_pitu)
pData(bg_pitu) <- data.frame(id=sampleNames(bg_pitu),group=c(0,0,0,1,1,1))
# 转录本水平的差异分析
stat_results_tran = stattest(bg,
feature='transcript',
meas='FPKM',getFC = TRUE,
covariate='group')
# 基因水平的差异分析
stat_results_gene_pitu = stattest(bg_pitu,
feature='gene',
meas='FPKM',getFC = TRUE,
covariate='group')
一些问题：
- 关于Deseq2处理之后出现的许多MSTRG.xxx的数据，我认为这些是新发现的转录本（非参考基因组对应的转录本），关于这一点推测的理由是：以羊为例参考基因组注释基因个数为34328，而经过Deseq2处理之后的数据会出现55353个"基因条目"，其中的有基因名字的有33509，另外的均为MSTRG.xxx数据，假如以不关注新的转录本为实验目的，可以将带基因命名的提取出来用于后续作图和分析。
- 关于Deseq2计算出的log2FC并非单纯的log2FC，其中包含可“一些缩小和缓和”的计算，但是大致可以近似为log2FC，比如(log2FC)=|1|代表FC的值的范围为FC>2&FC<0.5。
参考文章
- https://study.163.com/course/courseLearn.htm?courseId=1005231058#/learn/video?lessonId=1052094537&courseId=1005231058
- https://mp.weixin.qq.com/s/zNFvod8B-VhX7Kq7OgoRMA
- https://www.baidu.com/link?url=wxq34KJyLfCFq7k_L_4E6mJukZpe_s-ycjWM4X7Gu_TPosowriIuXNi1FFX4uw-lBsne5KbX434VrlH-3qxt2w4eiQtzktT5jbKkL05UTl3&wd=&eqid=cbf3fc9600059c0b000000056069657f
- https://www.bilibili.com/video/BV1KJ411p7WN
