一、序列比对
在2016年的一篇综述A survey of best practices for RNA-seq data analysis，提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。
1. Hisat2教程
1.1 下载安装
#conda直接安装
conda install hisat2
#源码下载安装
wget wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-source.zip
unzip hisat2-2.1.0-source.zip
make
1.2 构建index
直接下载现有的insex或通过Hisat2的方法进行创建
# 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites，但是先做的目的是为了提高运行效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index， 必须选项是基因组所在文件路径和输出的前缀
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19
1.3正式比对
hisat2基本用法就是hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> } [-S <hit>]，基本就是提供index的位置，PE数据或者是SE数据存放位置。然而其他可选参数却是进阶的一大名堂。新手就用默认参数呗。
hisat2 --dta -p 6 --max-intronlen 5000000 -x Oryza_sativa.IRGSP-1.0.genome -1 C1-1_good_1.fq -2 C1-1_good_2.fq -S C1-1.HISAT_aln.sam >hisat2_running.log 2>&1
1.4 Hisat2输出结果
比对之后会输出如下结果，解读一下就是全部数据都是100%的，2.88%的配对数据一次都没有比对，94.20%的数据比是唯一比对，2.92%是多个比对。然后如果不按照顺序来，有4.96%的比对。之后把剩下的部分用单端数据进行比对的话，65.57%数据没比对上，33.23%的数据比对一次，1.20%比对超过一次。零零总总的加起来是98.20%的比对。
20182824 reads; of these:
20182824 (100.00%) were paired; of these:
581893 (2.88%) aligned concordantly 0 times
19011569 (94.20%) aligned concordantly exactly 1 time
589362 (2.92%) aligned concordantly >1 times
----
581893 pairs aligned concordantly 0 times; of these:
28886 (4.96%) aligned discordantly 1 time
----
553007 pairs aligned 0 times concordantly or discordantly; of these:
1106014 mates make up the pairs; of these:
725197 (65.57%) aligned 0 times
367552 (33.23%) aligned exactly 1 time
13265 (1.20%) aligned >1 times
98.20% overall alignment rate
2. SAMtools三板斧
SAM（sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式，当然他可以用于存放未比对的数据。而目前处理SAM格式的工具主要是SAMTools
view: BAM-SAM/SAM-BAM 转换和提取部分比对
sort: 比对排序
merge: 聚合多个排序比对
index: 索引排序比对
faidx: 建立FASTA索引，提取部分序列
tview: 文本格式查看序列
#最常用的三板斧就是格式转换，排序，索引。
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
samtools index SRR35899${i}_sorted.bam
3. BAM/SAM文件格式
SAM文件主要由两个部分构成:
header：标记了该SAM文件的一些基本信息，比如版本、按照什么方式排序的、Reference信息等等。
本体：每行为一个reads，不同列记录了不同的信息，列与列之间通过tab分隔。
每列的含义：
MAPQ值：
表示为mapping的质量值，mapping Quality, It equals -10log10Pr{mapping position is wrong}, rounded to the nearest integer, A value 255 indicates that the mapping quality is not available. 该值的计算方法是mapping的错误率的-10log10值，之后四舍五入得到的整数，如果值为255表示mapping值是不可用的，如果是unmapped read则MAPQ为0，一般在使用bwa mem或bwa aln（bwa 0.7.12-r1039版本）生成的sam文件，第五列为60表示mapping率最高，一般结果是这一列的数值是从0到60，且0和60这两个数字出现次数最多
想把小于2的都丢弃：
samtools view -bSq 2 file.sam > filtered.bam
flag的含义：
1 ： 代表这个序列采用的是PE双端测序
2： 代表这个序列和参考序列完全匹配，没有插入缺失
4： 代表这个序列没有mapping到参考序列上
8： 代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
16：代表这个序列比对到参考序列的负链上
32 ：代表这个序列对应的另一端序列比对到参考序列的负链上
64 ： 代表这个序列是R1端序列， read1;
128 : 代表这个序列是R2端序列，read2；
256： 代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的
512： 代表这个序列在QC时失败了，被过滤不掉了（# 这个标签不常用）
1024: 代表这个序列是PCR重复序列（#这个标签不常用）
2048: 代表这个序列是补充的比对（#这个标签具体什么意思，没搞清楚，但是不常用）
cigar的含义:
cigar中会包含数字，代表了特定match持续了多少nt；以及不同的字符，代表了不同的match情况。
30S512M216N12S （30nt soft clip -> 512nt exact match -> 216nt skipped region -> 12nt soft clip）
30S (30nt soft clip)
40M (40nt exact match)
其中不同的字符及其含义如下：
参考：
https://www.jianshu.com/p/a584d31418f3
https://www.jianshu.com/p/9c87bba244d8
二、htseq-count的使用
HTSeq作为一款可以处理高通量数据的python包，由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人携手推出HTSeq — A Python framework to work with high-throughput sequencing data。自发布以来就备受广大分析人员青睐，其提供了许多功能给那些熟悉python的大佬们去自信修改使用，同时也兼顾着给小白们提供了两个可以拿来可用的可执行文件 htseq-count（计数） 和 htseq-qa（质量分析）。
具体参考：
https://www.cnblogs.com/triple-y/p/9338890.html
https://blog.csdn.net/herokoking/article/details/78257714
三、基因差异表达分析
1. DESeq2（DESeq2不支持无生物学重复的数据）
library("DESeq2")
#directory<-'E:\\R\\R-work'
directory<-'5.deseq2/15d_peiru'
sampleFiles<-grep('countf',list.files(directory),value=TRUE)
#sampleFiles
sampleCondition<-sub("(.*?)_.*.countf","\\1",sampleFiles)
sampleTable<-data.frame(sampleName=sampleFiles,fileName=sampleFiles,condition=sampleCondition)
#sampleTable
ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable,directory=directory,design=~condition)
#ddsHTSeq
colData(ddsHTSeq)$condition<-factor(colData(ddsHTSeq)$condition, levels=c("F1","za"))
dds<-DESeq(ddsHTSeq)
res<-results(dds)
res<-res[order(res$padj),]
#head(res)
png(file="plotMA.png",type="cairo")
plotMA(dds)
dev.off()
mcols(res,use.names=TRUE)
write.csv(as.data.frame(res),file='deseq2.csv')
png(file="plotDispEsts.png",type="cairo")
plotDispEsts(dds)
dev.off()
#pheatmap
sum(res$padj < 0.05, na.rm=TRUE)
library("pheatmap")
select <- order(rowMeans(counts(dds,normalized=TRUE)),decreasing=TRUE)[1:1000]
nt <- normTransform(dds) # defaults to log2(x+1)
log2.norm.counts <- assay(nt)[select,]
df <- as.data.frame(colData(dds))
pdf('heatmap.pdf',width = 6, height = 7)
pheatmap(log2.norm.counts, cluster_rows=TRUE, show_rownames=FALSE,cluster_cols=TRUE, annotation_col=df)
dev.off()
2. EdgeR
count.matrix
ID B.count wt.count
AT1G01010 384 314
AT1G01020 661 861
AT1G01030 48 47
AT1G01040 5608 6206
AT1G01046 77 82
AT1G01050 2889 2357
AT1G01060 6039 6296
AT1G01070 408 240
AT1G01073 0 0
edgeR.R
library(edgeR)
data = read.table('count.matrix', header=T, row.names=1, com='')
col_ordering = c(1,2)
rnaseqMatrix = data[,col_ordering]
rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(cpm(rnaseqMatrix) > 1) >= 2,]
conditions = factor(c(rep('B', 1), rep('wt', 1))) #对样本进行分组
exp_study = DGEList(counts=rnaseqMatrix, group=conditions) #建立基因表达列表，利用DEGList函数， 参数counts为read数矩阵，group为上一步的分组变量
exp_study = calcNormFactors(exp_study) #计算各样本内的标准化因子
exp_study_bcv = exp_study
bcv = 0.01 #估计生物学变异系数
et = exactTest(exp_study_bcv, dispersion = bcv ^ 2)
#gene = decideTestsDGE(et, p.value = 0.05, lfc = 0)
#summary(gene)
#head(res)
tTags = topTags(et,n=NULL)
result_table = tTags$table
result_table = data.frame(sampleA='B', sampleB='wt', result_table)
result_table$logFC = -1 * result_table$logFC
write.csv(result_table,file = 'B-wt.csv')
source('Rna-seq/trinityrnaseq-Trinity-v2.4.0/Analysis/DifferentialExpression/R/rnaseq_plot_funcs.R')
pdf('B-wt.pdf')
plot_MA_and_Volcano(result_table$logCPM, result_table$logFC, result_table$FDR)
dev.off()
网友评论