STEP4: 得到表达矩阵的流程

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了STEP4: 得到表达矩阵的流程相关的知识,希望对你有一定的参考价值。

参考技术A

这是RNA-Seq 上游分析的大致流程,比对+定量。当然实验目的若只需要定量已知基因,也可以选择free-alignment 的流程工具如kallisto/Salmon/Sailfish,其优点是可用于RNA-seq的基因表达的快速定量,但是对于小RNA和表达量低的基因分析效果并不好(2018年刚发表的一篇文章对free-alignment 的工具进行了质量评估,doi: https://doi.org/10.1101/246967 )。基于比对的流程,比对工具也有很多选择,如Hisat,STAR,Tophat(hista可以替代tophat),bowtie等, 还有据说速度超快的Subread。根据比对工具的选择,定量软件也可以配套选择,如
STAR/bowtie—RSEM;
Hisat —featureCounts/HTseq;
Subread—featureCounts。
由于后续要鉴定新的转录本,因此需要对转录本组装,组装可以选择cufflinks,或者stringtie。文章是用的Tophat+cufflinks。我的后续分析打算用Hisat2比对,stringtie组装,featureCounts定量,同时打算尝试下Subread+featureCounts的流程。
Hisat2+stringtie+featureCounts;
Subread+featureCounts

sra格式的数据需要先用fastq-dump转换, --split-3 表示双端测序,--gzip将生成的fastq文件压缩。

质控检测采用fastqc和multiqc联合使用, multiqc有利于多个样本同时展示质量检测结果。FastQC的检测报告左侧会给出测序质量的一个summary,通常并不是所有的检测项都会是绿色通过,会有一些警告和fail的项目。主要可以从Per base sequence quality 看一下测序碱基质量,Per sequence GC content 看一下GC含量,如果实际的GC含量(红线)出现双峰,且导致后期的序列比对很低时,需要引起注意,有可能是存在样品污染;再者就是看一下Adapter是否去除及去除的是否干净。这里的Adapter虽然显示通过,但是去除的并不干净,所以后续还会进一步去除Adapter。

Trimmomatic 支持多线程,处理数据速度快,主要用来去除 Illumina 平台的 fastq 序列中的接头,并根据碱基质量值对 fastq 进行修剪。软件有两种过滤模式,分别对应 SE 和 PE 测序数据,同时支持 gzip 和 bzip2 压缩文件。另外也支持 phred-33 和 phred-64 格式互相转化。

Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关,通常的过滤步骤如下:
ILLUMINACLIP : 过滤 reads 中的 Illumina 测序接头和引物序列,并决定是否去除反向互补的 R1/R2 中的 R2, 该引物序列可以在Trimmomatic软件的安装目录下找到,双端通常选择TruSeq3-PE-2。
SLIDINGWINDOW: 从 reads 的 5’ 端开始,进行滑窗质量过滤,切掉碱基质量平均值低于阈值的滑窗。
MAXINFO : 一个自动调整的过滤选项,在保证 reads 长度的情况下尽量降低测序错误率,最大化 reads 的使用价值。
LEADING : 从 reads 的开头切除质量值低于阈值的碱基。
TRAILING : 从 reads 的末尾开始切除质量值低于阈值的碱基。
CROP : 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
HEADCROP : 从 reads 的开头切掉指定数量的碱基。
MINLEN : 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
AVGQUAL : 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
TOPHRED33 : 将 reads 的碱基质量值体系转为 phred-33。
TOPHRED64 : 将 reads 的碱基质量值体系转为 phred-64。
参考 : NGS 数据过滤之 Trimmomatic 详细说明

Hisat2 比对,首先用hisat2-build 构建index,然后比对,输出sam格式,再用samtools将sam转为bam格式,并排序构建索引。
featureCounts 目前已经整合在subread,subread是一个综合的软件,还包括比对的工具和对应的R包, featrueCounts的定量效果和HTSeq差不多,但是featureCounts的优点非常快!

二进制版本的软件解压即可使用

(PS:存储不够了,后续选择两组数据做流程分析。)

SRR4042230和SRR4042231的比对率分别为88.02%和88.81%,总比对时间将近5小时。

hisat_counts.txt.summary 包含一些基本的统计信息。

可以看出subjunc比对不到2小时,比hisat2快3个多小时。

根据第1列是Geneid,第7,8列是counts数,用awk提取出geneID和counts。

一个植物转录组项目的实战
史上最快的转录组流程-subread
转录组定量-FEATURECOUNTS

Step By Step(Lua表达式和语句)

Step By Step(Lua表达式和语句)

一、表达式:

    1. 算术操作符:
    Lua支持常规算术操作符有:二元的“+”、“-”、“*”、“/”、“^”(指数)、“%”(取模),一元的“-”(负号)。所有这些操作符都可用于实数。然而需要特别说明的是取模操作符(%),Lua中对该操作符的定义为:
    a % b == a - floor(a / b) * b
    由此可以推演出x % 1的结果为x的小数部分,而x - x % 1的结果则为x的整数部分。类似的,x - x % 0.01则是x精确到小数点后两位的结果。
    
    2. 关系操作符:
    Lua支持的关系操作符有:>、<、>=、<=、==、~=,所有这些操作符的结果均为true或false。
    操作符==用于相等性测试,操作符~=用于不等性测试。这两个操作符可以应用于任意两个值。如果两个值的类型不同,Lua就认为他们不等。nil值与其自身相等。对于table、userdata和函数,Lua是通过引用进行比较的。也就是说,只有当他们引用同一个对象时,才视为相等。如:

 
1 a = {}
2 a.x = 1
3 a.y = 0
4 b = {}
5 b.x = 1
6 b.y = 1
7 c = a
 

    其结果是a == c,但a ~= b。
    对于字符串的比较,Lua是按照字符次序比较的。
   
    3. 逻辑操作符:
    Lua支持的逻辑操作符有:and、or和not。与条件控制语句一样,所有的逻辑操作符都将false和nil视为假,其他的结果均为真。和其他大多数语言一样,Lua中的and和or都使用“短路原则”。在Lua中有一种惯用写法"x = x or v",它等价于:if not x then x = v end。这里还有一种基于“短路原则”的惯用写法,如:
    max = (x > y) and x or y
    这等价于C语言中max = (x > y) ? x : y。由于x和y均为数值,因此它们的结果将始终为true。
    
    4. 字符串连接:
    前一篇Blog已经提到了字符串连接操作符(..),这里再给出一些简单的示例。
    /> lua
    > print("Hello " .. "World)
    Hello World
    > print(0 .. 1)  --即使连接操作符的操作数为数值类型,在执行时Lua仍会将其自动转换为字符串。
    01

    5. table构造器:
    构造器用于构建和初始化table的表达式。这是Lua特有的表达式,也是Lua中最有用、最通用的机制之一。其中最简单的构造器是空构造器{},用于创建空table。我们通过构造器还可以初始化数组,如:

 
 1 days = {"Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"}
2 for i = 1,#days do
3 print(days[i])
4 end
5 --输出结果为
6 --Sunday
7 --Monday
8 --Tuesday
9 --Wednesday
10 --Thursday
11 --Friday
12 --Saturday
 

    从输出结果可以看出,days在构造后会将自动初始化,其中days[1]被初始化为"Sunday",days[2]为"Monday",以此类推。
    Lua中还提供了另外一种特殊的语法用于初始化记录风格的table。如:a = { x = 10, y = 20 },其等价于:a = {}; a.x = 10; a.y = 20
    在实际编程时我们也可以将这两种初始化方式组合在一起使用,如:

 
polyline = {color = "blue", thickness = 2, npoints = 4, 
{x = 0, y = 0},
{x = 10, y = 0},
{x = -10, y = 1},
{x = 0, y = 1} }
print(polyline["color"]);
print(polyline[2].x)
print(polyline[4].y)
--输出结果如下:
--
blue
--
10
--
1
 

    除了以上两种构造初始化方式之外,Lua还提供另外一种更为通用的方式,如:

1 opnames = { ["+"] = "add", ["-"] = "sub", ["*"] = "mul", ["/"] = "div"}
2 print(opnames["+"])
3 i = 20; s = "-"
4 a = { [i + 0] = s, [i + 1] = s .. s, [i + 2] = s..s..s }
5 print(a[22])

    对于table的构造器,还有两个需要了解的语法规则,如:
    a = { [1] = "red", [2] = "green", [3] = "blue", }  
    这里需要注意最后一个元素的后面仍然可以保留逗号(,),这一点类似于C语言中的枚举。

    a = {x = 10, y = 45; "one", "two", "three" }
    可以看到上面的声明中同时存在逗号(,)和分号(;)两种元素分隔符,这种写法在Lua中是允许的。我们通常会将分号(;)用于分隔不同初始化类型的元素,如上例中分号之前的初始化方式为记录初始化方式,而后面则是数组初始化方式。

二、语句:

    1. 赋值语句:
    Lua中的赋值语句和其它编程语言基本相同,唯一的差别是Lua支持“多重赋值”,如:a, b = 10, 2 * x,其等价于a = 10; b = 2 * x。然而需要说明的是,Lua在赋值之前需要先计算等号右边的表达式,在每一个表达式都得到结果之后再进行赋值。因此,我们可以这样写变量交互:x,y = y,x。如果等号右侧的表达式数量少于左侧变量的数量,Lua会将左侧多出的变量的值置为nil,如果相反,Lua将忽略右侧多出的表达式。

    2. 局部变量与块:
    Lua中的局部变量定义语法为:local i = 1,其中local关键字表示该变量为局部变量。和全局变量不同的是,局部变量的作用范围仅限于其所在的程序块。Lua中的程序可以为控制结构的执行体、函数执行体或者是一个程序块,如:
    下面的x变量仅在while循环内有效。

1 while i <= x do
2 local x = i * 2
3 print(x)
4 i = i + 1
5 end

    如果是在交互模式下,当执行local x = 0之后,该变量x所在的程序即以结束,后面的Lua语句将被视为新的程序块。如果想避免此类问题,我们可以显式的声明程序块,这样即便是在交互模式下,局部变量仍然能保持其块内有效性,如:

1 do
2 local a2 = 2 * a
3 local d = (b ^ 2 - 4 * a) ^ (1 / 2)
4 x1 = (-b + d) / a2
5 x2 = (-b - d) / a2
6 end --a2和d的作用域至此结束。

    和其它编程语言一样,如果有可能尽量使用局部变量,以免造成全局环境的变量名污染。同时由于局部变量的有效期更短,这样垃圾收集器可以及时对其进行清理,从而得到更多的可用内存。    

    3. 控制结构:
    Lua中提供的控制语句和其它大多数开发语言所提供的基本相同,因此这里仅仅是进行简单的列举。然后再给出差异部分的详细介绍。如:
    1). if then else
    if a < 0 then
        b = 0
    else
        b = 1
    end
    
    2). if elseif else then
    if a < 0 then
        b = 0
    elseif a == 0 then
        b = 1
    else
        b = 2
    end
    
    3). while
    local i= 1
    while a[i] do
        print(a[i])
        i = i + 1
    end
    
    4). repeat
    repeat
        line = io.read()
    until line ~= "" --直到until的条件为真时结束。
    print(line)
    
    5). for
    for var = begin, end, step do --如果没有step变量,begin的缺省步长为1。
        i = i + 1
    end
    需要说明的是,for循环开始处的三个变量begin、end和step,如果它们使表达式的返回值,那么该表达式将仅执行一次。再有就是不要在for的循环体内修改变量var的值,否则会导致不可预知的结果。
    
    6). foreach
    for i, v in ipairs(a) do  --ipairs是Lua自带的系统函数,返回遍历数组的迭代器。
        print(v)
    end
    
    for k in pairs(t) do      --打印table t中的所有key。
        print(k)
    end
    见如下示例代码:

 
 1 days = {"Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday" }
2 revDays = {}
3 for k, v in ipairs(days) do
4 revDays[v] = k
5 end
6
7 for k in pairs(revDays) do
8 print(k .. " = " .. revDays[k])
9 end
10
11 --输出结果为:
12 --Saturday = 7
13 --Tuesday = 3
14 --Wednesday = 4
15 --Friday = 6
16 --Sunday = 1
17 --Thursday = 5
18 --Monday = 2
 

    7). break
    和C语言中的break语义完全相同,即跳出最内层循环。

以上是关于STEP4: 得到表达矩阵的流程的主要内容,如果未能解决你的问题,请参考以下文章

使用RSEM进行转录组测序的差异表达分析

R语言DESeq2基因差异表达分析

GSEA分析

Quantile Normalization

如何匹配,但排除正则表达式模式?

有向无环图——描述表达式