一、 首先放一段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的格式

一条序列有四行：

1. 序列综合信息
2. 序列碱基信息
3. +号
4. 序列对应的每个碱基的质量信息（根据荧光信号的强弱和弥散程度分析得到）
   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

主要用到的软件有：

1. 数据质量分析软件 fastqc multiqc
2. 数据质量控制软件 fastp(采用c++编写，功能强大，运行快)
3. 序列比对软件 hisat2（以及samtools将sam文件转化为bam文件）
4. 样本定量软件 stringtie
5. 差异表达分析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')

一些问题：

1. 关于Deseq2处理之后出现的许多MSTRG.xxx的数据，我认为这些是新发现的转录本（非参考基因组对应的转录本），关于这一点推测的理由是：以羊为例参考基因组注释基因个数为34328，而经过Deseq2处理之后的数据会出现55353个"基因条目"，其中的有基因名字的有33509，另外的均为MSTRG.xxx数据，假如以不关注新的转录本为实验目的，可以将带基因命名的提取出来用于后续作图和分析。
2. 关于Deseq2计算出的log2FC并非单纯的log2FC，其中包含可“一些缩小和缓和”的计算，但是大致可以近似为log2FC，比如(log2FC)=|1|代表FC的值的范围为FC>2&FC<0.5。

参考文章

1. https://study.163.com/course/courseLearn.htm?courseId=1005231058#/learn/video?lessonId=1052094537&courseId=1005231058
2. https://mp.weixin.qq.com/s/zNFvod8B-VhX7Kq7OgoRMA
3. https://www.baidu.com/link?url=wxq34KJyLfCFq7k_L_4E6mJukZpe_s-ycjWM4X7Gu_TPosowriIuXNi1FFX4uw-lBsne5KbX434VrlH-3qxt2w4eiQtzktT5jbKkL05UTl3&wd=&eqid=cbf3fc9600059c0b000000056069657f
4. https://www.bilibili.com/video/BV1KJ411p7WN