精华内容
下载资源
问答
  • RNA-seq数据分析实用方法,包含了RNA-seq数据分析的各个方面。
  • 本文基于文本数据与RNA-Seq数据在结构上具有的高度相似性,将文本数据分析中流行的概率模型LDA应用于RNA-Seq数据分析,设计了NU-LDA模型以测量读段在非均匀分布情况下基因和剪切异构体的表达水平。通过采用真实实验...
  • yeast_data:酵母RNA-seq数据分析
  • 半监督特征提取,用于RNA-Seq数据分析
  • 用于BioJupies项目( )的RNA-seq数据分析插件。 可在找到BioJupies Web服务器的源代码。 概述 什么是BioJupies? BioJupies是一个Web服务器,它使用户可以通过直观的界面从RNA-seq数据集自动生成Jupyter ...
  • RNA-Seq数据分析流程

    千次阅读 2020-11-03 13:09:43
    文章目录RNA-seq 数据分析流程相关软件安装下载数据sra转fastq格式数据质控数据质控,过滤低质量reads,去接头比对首先下载参考基因组及注释文件,建立索引比对sam文件转bam为bam文件建立索引reads的比对情况统计...

    RNA-seq 数据分析流程

    相关软件安装

    可以安装 conda,在后续其他软件安装时非常好用。可自行百度进行安装
    可根据文献调研,转录组数据分析所需软件列表:
    质控
    fastqc , multiqc, trimmomatic, cutadapt ,trim-galore
    比对
    star, hisat2, bowtie2, tophat, bwa, subread
    计数
    htseq, bedtools, deeptools, salmon

    在安装之前可以先 search一下安装包是否存在

    # conda search packagename
    conda install -y sra-tools
    conda install -y trimmomatic
    conda install -y cutadapt multiqc 
    conda install -y trim-galore
    conda install -y star hisat2 bowtie2
    conda install -y subread tophat htseq bedtools deeptools
    conda install -y salmon
    

    软件安装时,可以使用conda一步安装,也可以自行百度,下载源码,解压后添加环境变量,使用。
    但是在使用sra-tookit时出现了问题,因为已经提前用conda安装好了sra-tools,后面我在网页下载解压了sratoolkit,并且在环境变量中添加好了,导致两个软件版本不合,使用prefetch下载网页数据时总是报错。所以需要卸载其中一种,保持版本一致。
    自行安装sratoolkit可参考:sratoolkit安装

    下载数据

    进入网页,得到SRR_Acc_list.txt

    SRR957677
    SRR957678
    SRR957679
    SRR957680
    

    批量下载

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    cat SRR_Acc_List.txt | while read id; do (nohup prefetch  ${id} &);done
    #ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done #在内地下载速度很慢,所以我杀掉这些下载进程
    

    当显示所有下载都sucessful时,可以进行下一步
    在这里插入图片描述
    下载完成之后,文件会自动保存在目录:~/ncbi/public/sra/

    mv ~/ncbi/public/sra/*.sra $wkd/rawdata
    

    sra转fastq格式

    使用fastq-dump

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    
    for i in $wkd/rawdata/*sra
    do
            echo $i
            fastq-dump --split-3 --skip-technical --clip --gzip $i  ## 批量转换
    done
    

    原始数据如果是双端测序结果,fastq-dump配合–split-3参数,一个样本被拆分成两个fastq文件,如果是单端测序,就只生成一个文件

    数据质控

    fastqc生成质控报告,multiqc将各个样本的质控报告整合为一个。

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    
    ls $wkd/rawdata/*gz | xargs fastqc -t 2
    multiqc ./ 
    

    得到结果如下:

    ├── [4.0K]  multiqc_data
    │   ├── [2.1M]  multiqc_data.json
    │   ├── [6.8K]  multiqc_fastqc.txt
    │   ├── [2.2K]  multiqc_general_stats.txt
    │   ├── [ 16K]  multiqc.log
    │   └── [3.4K]  multiqc_sources.txt
    ├── [1.5M]  multiqc_report.html
    ├── [236K]  SRR1039508_1_fastqc.html
    ├── [279K]  SRR1039508_1_fastqc.zip
    ├── [238K]  SRR1039508_2_fastqc.html
    ├── [286K]  SRR1039508_2_fastqc.zip
    

    每个id_fastqc.html都是一个质量报告,multiqc_report.html是所有样本的整合报告

    数据质控,过滤低质量reads,去接头

    如果双端测序,运行一下代码,得到两列数据,config文件

    mkdir $wkd/clean 
    cd $wkd/clean 
    ls /home/jmzeng/project/airway/raw/*_1.fastq.gz >1
    ls /home/jmzeng/project/airway/raw/*_2.fastq.gz >2
    paste 1 2  > config
    

    并打开qc.sh,写入一下内容:

    bin_trim_galore=trim_galore
    dir='/home/jmzeng/project/airway/clean'
    cat $1 |while read id
    do
            arr=(${id})
            fq1=${arr[0]}
            fq2=${arr[1]} 
     $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir  $fq1 $fq2 
    done 
    

    单端测序:

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    mkdir $wkd/cleandata 
    cd $wkd/cleandata 
    ls /home/meiling/baiduyundisk/RNA-seq/rawdata/*.fastq.gz >config
    #ls /home/jmzeng/project/airway/raw/*_2.fastq.gz >2
    #paste 1 2  > config
    

    qc.sh

    wkd=/home/meiling/baiduyundisk/RNA-seq/cleandata #设置工作目录
    cd $wkd
    cat config |while read id
    do
            arr=(${id})
            fq1=${arr[0]}
    #       fq2=${arr[1]} 
        trim_galore -q 25 --phred33 --length 36 --stringency 3 -o $wkd  $fq1 
    done
    

    其他软件,对于单端测序和双端测序的不同用法,可以自行百度具体软件用法,都有解释,参考官方参数即可。

    比对

    首先下载参考基因组及注释文件,建立索引

    基因组相关内容大家可以自行百度,但是参考基因组和注释文件一定要对应在同一个网站,同一个版本。可以下载的网站有UCSC,NCBI,ENSEMBEL
    我是在ENSEMBEL上面下载的,一般选择primary参考基因组,注释文件选择scaff
    复制链接地址后,axel 下载,非常快。下载后解压。

    axel -n 20 ftp://ftp.ensembl.org/pub/release-101/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
    axel -n 20 ftp://ftp.ensembl.org/pub/release-101/gtf/homo_sapiens/Homo_sapiens.GRCh38.101.chr_patch_hapl_scaff.gtf.gz
    ls
    #Homo_sapiens.GRCh38.101.chr_patch_hapl_scaff.gtf.gz
    #Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
    ##解压
    gzip -d Homo_sapiens.GRCh38.101.chr_patch_hapl_scaff.gtf.gz
    gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
    ls -lh
    total 4.3G
    drwxr-xr-x 2 meiling meiling 4.0K Nov  2 21:26 hisat2_index
    -rw-r--r-- 1 meiling meiling 1.3G Nov  1 20:48 Homo_sapiens.GRCh38.101.chr_patch_hapl_scaff.gtf
    -rw-r--r-- 1 meiling meiling 3.0G Nov  1 20:50 Homo_sapiens.GRCh38.dna.primary_assembly.fa
    ##建立索引
    hisat2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa hisat2_index/hisat_hg38
    

    在这里插入图片描述

    比对

    star, hisat2, bowtie2, tophat, bwa, subread都是可以用于比到的软件
    批量比对:

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    
    cd $wkd/cleandata 
    ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
    ##单端测序
    ls -lh ${id}_trimmed.fq.gz
    hisat2 -q -p 10 -x $wkd/reference/hg38/hisat2_index/hisat_hg38 -U ${id}_trimmed.fq.gz -S $wkd/align/${id}.hisat.sam
    ##双端测序
    #ls -lh ${id}_1_val_1.fq.gz   ${id}_2_val_2.fq.gz 
    #hisat2 -p 10 -x $wkd/reference/hg38/hisat2_index/hisat_hg38 -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -S ${id}.hisat.sam
    # subjunc -T 5  -i /public/reference/index/subread/hg38 -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fq.gz -o ${id}.subjunc.sam
    # bowtie2 -p 10 -x /public/reference/index/bowtie/hg38  -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -S ${id}.bowtie.sam
    # bwa mem -t 5 -M  /public/reference/index/bwa/hg38   ${id}_1_val_1.fq.gz   ${id}_2_val_2.fq.gz > ${id}.bwa.sam
    done
    

    这里是演示多个比对工具,但事实上,对RNA-seq数据来说,不要使用bwa和bowtie这样的软件,它需要的是能进行跨越内含子比对的工具。
    比对结果输出:

    bash align.sh > align.log 2>&1
    ##打开align.log
    ###单端测序比对结果输出
    -rw-r--r-- 1 meiling meiling 984M Nov  2 16:33 SRR957677_trimmed.fq.gz
      File "/opt/software/hisat2-2.2.0/hisat2_read_statistics.py", line 182
        length_map = sorted(length_map.iteritems(), key=lambda (k,v):(v,k), reverse=True)
                                                               ^
    SyntaxError: invalid syntax
    20498068 reads; of these:
      20498068 (100.00%) were unpaired; of these:
        689901 (3.37%) aligned 0 times
        17438548 (85.07%) aligned exactly 1 time
        2369619 (11.56%) aligned >1 times
    96.63% overall alignment rate
    

    在这里插入图片描述

    sam文件转bam

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    cd $wkd/align 
    ls *.sam|cut -d"." -f 1 |while read id ;do
        samtools view -@ 8 -S {$id}.sam -1b -o {$id}.bam
        samtools sort -@ 8 -l 5 -o {$id}.sort.bam {$id}.bam
    done
    #rm *.sam 
    

    在这里插入图片描述

    为bam文件建立索引

    ls *.sort.bam |xargs -i samtools index {}
    

    reads的比对情况统计

    ls *.sort.bam |while read id ;do ( samtools flagstat -@ 1 $id >  $(basename ${id} ".bam").flagstat  );done
    

    在这里插入图片描述

    计数 counts

    首先安装featurecounts,这个软件已经整合到subread里面了

    axel -n 20 https://udomain.dl.sourceforge.net/project/subread/subread-2.0.1/subread-2.0.1-Linux-x86_64.tar.gz
    tar -zxvf subread-2.0.1-Linux-x86_64.tar.gz 
    cd subread-2.0.1-Linux-x86_64/
    cd bin
    pwd
    #添加环境变量
    

    计数

    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    cd $wkd/align 
    # 如果一个个样本单独计数,输出多个文件使用代码是:
    for fn in {508..523}
    do
    featureCounts -T 5 -p -t exon -g gene_id  -a $gtf -o $fn.counts.txt SRR1039$fn.bam
    done
    # 如果是批量样本的bam进行计数,使用代码是:
    wkd=/home/meiling/baiduyundisk/RNA-seq #设置工作目录
    cd $wkd/align 
    # 如果一个个样本单独计数,输出多个文件使用代码是:
    # for fn in {508..523}
    # do
    # featureCounts -T 5 -p -t exon -g gene_id  -a /public/reference/gtf/gencode/gencode.v25.annotation.gtf.gz -o $fn.counts.txt SRR1039$fn.bam
    # done
    # 运行featureCounts对所有样本进行基因水平定量
    # 如果是批量样本的bam进行计数,使用代码是:
    # 设置featureCounts所需要的gtf文件位置(注释文件)
    gtf="$wkd/reference/hg38/Homo_sapiens.GRCh38.101.chr_patch_hapl_scaff.gtf.gz"   
    featureCounts -T 5 -t exon -g gene_id  -a $gtf -o  $wkd/counts/all.id.txt  *.sort.bam  1>counts.id.log 2>&1
    #解释:-T 1:线程数为1;-p:表示数据为paired-end,双末端测序数据;-t exon表示 feature名称为exon;-g gene_id表示meta-feature名称为gene_id(ensembl名称);-a $gtf 表示输入的GTF基因组注释文件;-o all.id.txt 设置输出文件类型;*.sort.bam是被分析的bam文件。featureCounts支持通配符*
    # 这样得到的  all.id.txt  文件就是表达矩阵啦,但是,这个 featureCounts有非常多的参数可以调整。
    cd $wkd/counts 
    # 使用multiqc对featureCounts统计结果进行可视化。
    multiqc all.id.txt.summary
    # 只保留all.id.txt文件的【基因名】和【样本counts】
    cat all.id.txt | cut -f 1,7 |grep -v '^#' > counts.txt
    

    运行结束:
    在这里插入图片描述
    得到以下文件:
    在这里插入图片描述

    # Program:featureCounts v2.0.1; Command:"featureCounts" "-T" "5" "-t" "exon" "-g" "gene_id" "-a" "/home/meiling/baiduyundisk/RNA-seq/reference/hg38/Homo_sapiens.GRCh38.101.chr_patch_hapl_scaff.gtf.gz" "-o" "/home/meiling/baiduyundisk/RNA-seq/counts/all.id.txt" "SRR957677.sort.bam" "SRR957678.sort.bam" "SRR957679.sort.bam" "SRR957680.sort.bam" 
    Geneid	Chr	Start	End	Strand	Length	SRR957677.sort.bam	SRR957678.sort.bam	SRR957679.sort.bam	SRR957680.sort.bam
    ENSG00000223972	1;1;1;1;1;1;1;1;1	11869;12010;12179;12613;12613;12975;13221;13221;13453	12227;12057;12227;12721;12697;13052;13374;14409;13670	+;+;+;+;+;+;+;+;+	1735	0	0	0	0
    ENSG00000227232	1;1;1;1;1;1;1;1;1;1;1	14404;15005;15796;16607;16858;17233;17606;17915;18268;24738;29534	14501;15038;15947;16765;17055;17368;17742;18061;18366;24891;29570	-;-;-;-;-;-;-;-;-;-;-	1351	7	1	2	3
    ENSG00000278267	1	17369	17436	-	68	0	0	0	0
    ENSG00000243485	1;1;1;1;1	29554;30267;30564;30976;30976	30039;30667;30667;31109;31097	+;+;+;+;+	1021	0	0	0	0
    ENSG00000284332	1	30366	30503	+	138	0	0	0	0
    ENSG00000237613	1;1;1;1;1	34554;35245;35277;35721;35721	35174;35481;35481;36081;36073	-;-;-;-;-	1219	0	0	0	0
    ENSG00000268020	1	52473	53312	+	840	0	0	0	0
    ENSG00000240361	1;1;1;1	57598;58700;62916;62949	57653;58856;64116;63887	+;+;+;+1414	0	0	0	0
    

    可以单个样本计数,之后在进行合并,使用Python脚本或者R语言merge
    具体可以参考这里
    或者直接批量比对,输出的all.id.txt即为合并后的表达矩阵。可以在这里选择注释,或者找完差异基因后在进行注释。

    差异基因分析

    后续需要在R里面进行差异基因分析,DESeq和edger里面自带标准化的函数,后续还可以进行KEGG,GSEA,GO富集分析

    展开全文
  • RNAdetector:免费的用户友好型独立软件,用于RNA-Seq数据分析 RNAdetector是一种用户友好的软件,可根据RNA-Seq数据分析蛋白质编码基因和ncRNA。 进行分析 RNAdetector可以执行多种类型的分析,例如: 量化和归一...
  • BingleSeq-用于批量和单细胞RNA-Seq数据分析的用户友好型R包。 手稿可作为。 安装 可以使用以下代码直接从GitHub安装BingleSeq: library("devtools") install_github("dbdimitrov/BingleSeq") # Start the ...
  • RNA-seq数据分析(HISAT2+featureCounts+StringTie)

    千次阅读 多人点赞 2020-08-21 21:53:58
    RNA-seq数据分析简介 简介 基因表达是功能基因组学研究的一个重要领域。基因表达与基因信息从基因组DNA模板到功能蛋白产物的流动有关(图1)。大规模并行RNA测序(RNA-seq)已成为一种标准的基因表达检测方法,尤其用于...


    简介

    基因表达是功能基因组学研究的一个重要领域。基因表达与基因信息从基因组DNA到功能蛋白产物的流动有关。RNA-seq已成为一种标准的基因表达检测方法,尤其用于检测转录本相对丰度和多样性。已有研究表明,其可靠性可以与其他成熟的方法如定量聚合酶链式反应(qPCR)相媲美。
    声明: 文中如有不足请留言讨论!

    1 生物基础

    吐槽: 系统性总结真的好困难,刚开始阅读文献也很吃力。但是行动是知识特有的果实,慢慢积累吧。

    1.1 中心法则

    中心法则对于学习生物的人来说再熟悉不过了,基因信息的流动,蛋白的产生等一系列生物学事件都基于中心法则,是基因表达分析的一条主线。遗传信息从双链基因组DNA模板到翻译后修饰蛋白的流动是每个阶段至关重要的分子特征。RNA-seq通常针对成熟的mRNA分子。

    在这里插入图片描述
    图中对于中心法则做了系统的总结,感兴趣可以click here

    1.2 RNA-seq Protocol

    典型的RNA-seq流程包括从感兴趣的样本中分离RNA、建库、高通量测序产生数以百万计的reads(长度一般为30-300bp)、比对到参考基因组或转录组,差异表达分析、发现转录本亚型和其它的下游分析。下图很好的概括了RNA-seq数据产生的过程。

    在这里插入图片描述
    RNA-seq的应用:
    检测所有基因在特定条件下的表达水平(发育阶段、不同组织、正常与疾病、药物治疗)。 发现新基因,可变剪切,基因突变和基因融合等。

    1.3 RNA-seq总的路线图

    这张图引自A survey of best practices for RNA-seq data analysis。包括了三个板块,分别是预分析,核心分析和高级分析。很好的概括了实验设计和拿到测序reads后的数据分析,并介绍了不同的分析路线和每一步的数据分析方法,可以说是一个很好的大纲。
    在这里插入图片描述
    更多的细节可以通过阅读文章中的引用进一步了解。

    2 数据分析

    注意:

    1. 这里省略了进入相关文件目录的步骤,大家分析时要注意。
    2. 本次数据是小鼠早期胚胎测序得到的双末端数据,单末端相关参数可以使用–help查看帮助文档。但分析总体流程相同。
      我使用的是HISAT2+featureCounts+StringTie流程。

    2.1 前期准备

    2.1.1 创建目录&安装conda

    创建目录

    # 首先使用cd命令需要进入自己的目录
    cd ~/
    # 创建文件夹,用于存储参考基因组
    mkdir -p /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/
    # 建立软件安装目录
    mkdir biosoft
    # 此外还需要建立项目分析的目录以及备份文件,方便查找
    

    conda安装

    Conda 是一种通用包管理系统,旨在构建和管理任何语言的任何类型的软件。通常与 Anaconda (集成了更多软件包,https://www.anaconda.com/download/#download) 和 Miniconda(只包含基本功能软件包, https://conda.io/miniconda.
    html) 一起分发。

    # 从官网下载最新版Miniconda3安装包,但速度较慢
    wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh 
    # 安装Miniconda3:安装过程只需要输入 yes 或者按 Enter
    bash Miniconda3-latest-Linux-x86_64.sh
    # miniconda3安装成功,并成功配置好环境变量
    source ~/.bashrc
    # 镜像设置
    # 下面这四行配置清华大学的bioconda的channel地址,国内用户推荐
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
    # 配置镜像成功
    conda config --set show_channel_urls yes
    conda config --set channel_priority flexible
    # 查看配置文件
    cat ~/.condarc
    

    在这里插入图片描述

    conda常用命令

    conda create -y -n trans python=3 # 创建小环境并成功安装python3
    conda info --e # 查看当前conda环境
    conda-env list # 列出所有小环境
    conda activate trans # 激活小环境
    conda deactivate # 退出小环境
    

    2.1.2 常用文件格式简介

    1. SRA:NCBI SRA数据库存放格式
    SRA(Sequence Read Archive):SRA是一个数据库,NCBI为了解决高通量数据庞大的存储压力而设计的一种数据压缩方案。
    一般使用fastq-dumpfasterq-dump来将其转换为fastq格式的数据,才能做后续分析。
    2. FASTQ:高通量数据存储格式
    FASTQ格式将序列和Phred质量存储在一个文件中。序列和质量得分皆由单个ASCII字符编码。最早由Sanger Institute开发使用,目前已经是高通量测序结果的标准格式。
    在这里插入图片描述
    在FASTQ文件中,一个序列通常由四行组成:
    第一行由@开头,之后为序列的标识符和描述信息;
    第二行为序列信息;
    第三行以+开头,之后可以再次加上序列的标识和描述信息;
    第四行为质量的分信息,与第二行的序列相对应,其长度必须与第二行相同。

    Phred质量值简介:
    在这里插入图片描述
    3. FASTA:记录信息的格式
    对于每条序列,首行以“>”开头,其之后是注释;
    在首行之后,是以单字母标准编码表达的实际序列数据
    FASTA的序列表达:
    1) 核算编码:A,C,G,T,N,U,R,Y,K,M,S,W,B,D,H,V
    2) 氨基酸编码:[A-Z],*(代表终止密码子)

    fasta序列格式如下:
    在这里插入图片描述
    4. SAM/BAM:高通量数据比对存放格式
    SAM文件是由比对产生的以tab建分割的文件格式,BAM是SAM文件的二进制压缩版本。使用samtools view -S -b -o my.bam my.sam可以将SAM文件转换为BAM文件。
    在这里插入图片描述
    5. BED:基因组浏览器常用格式
    BED 文件格式提供了一种灵活的方式来定义的数据行,以用来描述注释信息。BED行有3个必须的列和9个额外可选的列。每行的数据格式要求一致。

    必须包含的3列:

    (1). chrom - 染色体名字(e.g. chr3,chrY, chr2_random)或scafflold 的名字(e.g. scaffold0671 )。
    (2). chromStart - 染色体或scaffold的起始位置,染色体第一个碱基的位置是0。
    (3). chromEnd - 染色体或scaffold的结束位置,染色体的末端位置没有包含到显示信息里面。例如,首先得100个碱基的染色体定义为chromStart =0 .chromEnd=100, 碱基的数目是0-99。

    9 个额外的可选列:
    (4). name - 指定BED行的名字,这个名字标签会展示在基因组浏览器中的bed行的左侧。
    (5). score - 0到1000的分值,如果在注释数据的设定中将原始基线设置为1,那么这个分值会决定现示灰度水平(数字越大,灰度越高),下面的这个表格显示GenomeBrowser

    shade	 	 	 	 	 	 	 	 	score in range  	≤ 166	167-277	278-388	389-499	500-611	612-722	723-833	834-944	≥ 945
    (6). strand - 定义链的方向,’’+” 或者”-”。
    (7). thickStart - 起始位置(The starting position atwhich the feature is drawn thickly)(例如,基因起始编码位置)。
    (8). thickEnd - 终止位置(The ending position at whichthe feature is drawn thickly)(例如:基因终止编码位置)。
    (9). itemRGB - 是一个RGB值的形式, R, G, B (eg. 255, 0,0), 如果itemRgb设置为’On”, 这个RBG值将决定数据的显示的颜色。
    (10). blockCount - BED行中的block数目,也就是外显子数目。
    (11). blockSize - 用逗号分割的外显子的大小, 这个item的数目对应于BlockCount的数目。
    (12). blockStarts - 用逗号分割的列表, 所有外显子的起始位置,数目也与blockCount数目对应。
    6. GFF/GTF
    GFF(General Feature Format)是基于Sanger GFF2的一种格式。GFF有9个必需字段,这些字段必须用制表符分隔。如果字段用空格而不是制表符分隔,则将不能正确显示。GTF (Gene Transfer Format, GTF2.2)是GFF的一种扩展格式。
    在这里插入图片描述
    在这里插入图片描述

    2.2 软件安装

    这里介绍三种软件安装方法。

    2.2.1 conda安装软件

    # 可以一次安装多个软件
    conda install -y -c hcc aspera-cli
    conda install -y sra-tools
    conda install -y fastqc trim-galore hisat2 multiqc samtools featureCounts
    # 运行以下命令,检查软件是否安装成功
    hisat2 --help
    which ascp # 查看软件安装路径
    

    2.2.2 预编译版本软件安装

    # 以aspera为例
    wget -c https://d3gcli72yxqn2z.cloudfront.net/connect_3_9_1_171801_ga/bin/ibm-aspera-connect-3.9.1.171801-linux-g2.12-64.tar.gz
    # 解压并运行脚本
    tar -xzf ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz  
    bash ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.sh
    # 添加环境变量
    export PATH=~/.aspera/connect/bin:\$PATH >>~/.bashrc 
    # source使变量生效
    source ~/.bashrc
    

    2.2.3 源码安装

    # 首先使用wget下载
    # 安装分三步:
    # 配置 
    ./configure --prefix=想要指定的路径
    # 编译
    make
    # 安装
    make install
    # 具体的安装在下载之后可以查看readme文件,一般有详细的介绍
    

    2.3 数据下载

    我这里使用aspera下载sra数据。
    首先我们通过阅读文章获取数据的GEO accession,在GEO数据库中找到对应的BioProject,这里建议在ebi数据库下载对应的SRR号的相关信号,并获取aspera下载链接以及md5值等相关信息。

    # ------------------------------------------------------
    # GEO accession:GSE70605
    # ------------------------------------------------------
    # 下载得到tsv文件首先需要使用awk命令提取文件内aspera链接,保存为srr.aspra文件。之后运行下面脚本
    创建一个脚本,命名为aspera.sh:
    
    ```bash
    # 找到openssh文件位置,替换以下代码文件路径
    find /home/hwb/ -name asperaweb_id_dsa.openssh
    
    #!/bin/bash
    cat srr.aspera | while read id
    do 
    ascp -QT -l 300m -P33001 -i ~/miniconda3/envs/aspera/etc/asperaweb_id_dsa.openssh era-fasp@${id} ./
    done
    
    # 运行脚本:
    ./aspera.sh
    

    下载完成之后一定要检查数据完整性,不然分析白做

    使用awk命令提取md5值,生成md5文件,运行以下命令检查文件完整性

    md5sum -c md5

    md5文件如下图:
    ![在这里插入图片描述](https://img-blog.csdnimg.cn/20200821200722787.PNG#pic_center)
    ## 2.4 sra2fastq
    将sra文件转换为fastq文件
    ```bash
    创建脚本文件,个人觉得这样很方便。方便下次使用
    #!bin/bash
    cat srr.list | while read id
    do 
    nohup fastq-dump --split-3 $id -O ./ &
    done
    

    2.5 质控(QC)

    这一步生成的是HTML文件,可以用浏览器打开查看。其中,各种颜色是各项标准分析结果:绿色代表"PASS";黄色代表"WARN";红色代表"FAIL"。

    ls ../raw/*.fastq|xargs fastqc -t 10 -o ./
    multiqc ./# 整合质控结果 
    

    在这里插入图片描述

    2.6 去接头

    去接头完成之后,一定要再次质控,确保数据tidy

    同样,写个简单的脚本
    #!/bin/bish
    # paired sequence
    # trim adapter
    paste <(ls *1.fastq) <(ls *2.fastq) > config
    # mkdir ../clean
    cat config | while read id
    do
    arr=$id
    fq1=${arr[0]}
    fq2=${arr[1]}
    nohup trim_galore --gzip --paired $fq1 $fq2 -o ../clean > ../clean/trim.log &
    done
    

    2.7 hisat2比对

    我这里使用的是HISAT2,index是在HISAT2官网下载的。所以省去了建index的步骤。hisat2的index有八个文件。这里建议新手将参考基因组的fa文件,注释gtf文件和index文件放在同一文件夹下,不然报错很难找。
    在这里插入图片描述

    这里也是一个小脚本
    #!/bin/bash
    # hisat2 of paired aligned
    # dir = pwd
    # 制作文件,第一列为输出文件名,第二列为双末端数据第一个文件,第三列为第二个文件
    paste <(ls *fq.gz | cut -d"_" -f 1 |sort | uniq) <(ls *1_val_1.fq.gz) <(ls *2_val_2.fq.gz) > configure
    # mkdir ../sam_out
    cat configure|while read id;do
    arr=(${id})
    sample=${arr[0]}
    fq1=${arr[1]}
    fq2=${arr[2]}
    echo $sample $fq1 $fq2
    nohup hisat2 -p 4 -k 1 --dta -x /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/genome_tran -1 $fq1 -2 $fq2 -S ../sam
    _out/$sample.sam &
    done
    

    到此生成了sam文件,比对的结果文件,这个文件一般较大,我们需要将其转换为bam文件,以便下一步分析。当然,可以在这个脚本中再添加命令,直接生成bam文件。

    # sam2bam
    cd ../sam_out/
    ls *sam|while read id;do nohup samtools view -b -F 4 -@ 4 $id | samtools sort -@ 4 -O BAM - > ../bam_out/$(basename $id "sam")mapped.sort.bam &
    done
    

    2.8 定量和转录本组装

    raw counts:
    featuresCounts软件用于统计基因/转录本上mapping的reads数,也就是用于raw count定量。该软件不仅支持基因/转录本的定量,也支持exon, gene bodies, genomic bins, chromsomal locations等区间的定量。此外,还有htseq-count 等。sam和bam文件均可作为输入文件。

    # 单样本定量
    featureCounts
    -T 5  \
    -t exon \
    -g gene_id -a /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/Mus_musculus.GRCm38.84.chr.gtf -o counts.txt mapping.sam
    # 多样本
    featureCounts -t exon -g gene_id -a /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/Mus_musculus.GRCm38.84.chr.gtf -o counts.txt ./*.bam
    

    raw count矩阵文件简单处理,就可以用R包DESeq2筛选差异基因,可以参考我的另一篇博客,转录组差异表达分析和火山图可视化
    StrintTie进行转录本组装和定量:
    gtf结果文件中有coverage,TPM和FPKM。此外RSEM也可计算FPKM值。

    cd ../bam
    ls *bam|while read id;do stringtie $id -p 4 -G /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/Mus_musculus.GRCm38.84.chr.gtf -o ../out_gtf/$(basename $id "bam")gtf;done
    

    至此,RNA-seq上游数据分析基本已经完成。当然,有许多的分析流程,比如salmon流程,tophat+cufflinks等。使用之前尽量比对各个流程的优缺点。

    数据分析过程坑很多,踩过了才能真正入门啊!(个人观点)


    最后,非常感谢Jimmy老师的视频分享! 由于篇幅太长,本文对于软件的使用没有详细介绍,大家可以去相应软件官网或者使用–help命令详细了解。


    ##如有侵权请联系作者删除!

    欢迎关注我的微信公众号。在这里插入图片描述

    参考

    [1] 生信技能树jimmy老师一系列课程

    [2] A survey of best practices for RNA-seq data analysis

    [3] Informatics for RNA Sequencing: A Web Resource for Analysis on the Cloud

    展开全文
  • 开放的RNA-Seq数据分析管道教程,其中包含最近处理Zika病毒研究中的数据的示例 辰王和AVI Ma'ayan BD2K-LINCS数据协调和集成中心(DCIC)位于纽约西奈山的伊坎医学院,纽约州10029美国 出版物 Wang Z和Ma'ayanA。一...
  • 生信入门(三)——使用DESeq2进行RNA-seq数据分析 文章目录生信入门(三)——使用DESeq2进行RNA-seq数据分析一、学习目标二、实验数据1、数据来源2、建模计数数据3、转录本丰度4、salmon定量三、用tximport导入R1...

    生信入门(四)——使用DESeq2进行RNA-seq数据分析


    今日学习内容DESeq2分析RNA-seq数据

    一、学习目标

    • 直观地评估 RNA-seq 数据的质量
    • 对 RNA-seq 数据进行基本的差异分析
    • 将 RNA-seq 结果与其他实验数据进行比较
    • 从 FASTQ 文件中量化转录表达
    • 将定量导入 R/Bioconductor
    • 执行质量控制和探索性数据分析
    • 执行差异表达
    • 与其他实验数据重叠
    • 构建动态报告
    • 整合 DESeq2 和 zinbwave 以获取单细胞 RNA-seq 数据

    二、实验数据

    1、数据来源

    GEO 条目 GSE52778

    • 本文将展示如何使用替代数据集(tximport小插图中使用的tximportData包)导入 RNA-seq 量化数据。后来我们将加载了数道数据集,其使用计数summarizeOverlaps从GenomicAlignments包。如下所述,我们推荐使用tximport管道来生成计数矩阵,但我们还没有包含气道数据集必要量化文件的 Bioconductor 包。

    2、建模计数数据

    作为输入,基于计数的统计方法,例如DESeq2 , edgeR , limma with the voom method , DSS、EBSeq 和baySeq期望输入数据是从 RNA-seq 或其他高通量测序实验中获得的计数矩阵。第i行和j 中的值矩阵的第 -th 列表示样本j 中的基因i分配了多少读数(或片段,对于双端 RNA-seq)。类似地,对于其他类型的检测,矩阵的行可能对应于例如结合区域(使用 ChIP-Seq)、细菌种类(使用宏基因组数据集)或肽序列(使用定量质谱)。

    3、转录本丰度

    在此工作流程中,我们将展示如何使用由Salmon (Patro et al. 2017)软件包量化的转录本丰度。Salmon和其他方法,例如Sailfish 、kallisto 或RSEM ,估计所有(已知的、带注释的)转录本的相对丰度,而没有对齐读取。因为估计转录本的丰度涉及推理步骤,所以估计计数. 大多数方法要么使用称为估计最大化或贝叶斯技术的统计框架来估计丰度和计数。在量化之后,我们将使用tximport (Soneson、Love 和 Robinson 2015)包来组装估计计数和偏移矩阵,以便与 Bioconductor 差异基因表达包一起使用。

    将转录本丰度量词与tximport结合使用来生成基因级计数矩阵和归一化偏移的优点是:

    • 这种方法纠正了样本中基因长度的任何潜在变化(例如,来自不同异构体的使用)
    • 与基于对齐的方法相比,其中一些方法要快得多,需要更少的内存和磁盘使用量
    • 可以避免丢弃那些可以与具有同源序列的多个基因对齐的片段(Robert and Watson 2015)。

    请注意,转录本丰度量词会跳过存储读取比对(SAM 或 BAM 文件)的大文件的生成,而是生成存储每个转录本的估计丰度、计数和有效长度的较小文件。
    salmon量化转录本丰度教程

    4、salmon定量

    略~因为我看不懂
    在这里插入图片描述

    三、用tximport导入R

    1、指定文件位置

    量化之后,我们可以使用tximport将数据导入 R 并使用 Bioconductor 包进行统计分析。通常,我们只需将tximport指向quant.sf我们机器上的文件。但是,因为我们将这些文件作为 R 包的一部分进行分发,所以我们必须执行一些额外的步骤,以确定 R 包以及这些文件在您的计算机上的位置。
    R 函数system.file可用于找出包中的文件安装在计算机上的哪个位置。这里我们要求提供R 包存储外部数据的目录的完整路径,这是tximportData包的一部分。

    BiocManager::install("tximportData")
    library("tximportData")
    dir<-system.file("extdata",package = "tximportData")
    list.files(dir)
    list.files(file.path(dir,"salmon"))#量化后的数据在salmon中
    

    在这里插入图片描述

    • 此处使用的标识符是来自European Nucleotide Archive的ERR标识符。

    我们需要创建一个指向量化文件的命名向量。我们将首先通过读取包含样本 ID 的表来创建文件名向量,然后将其与和结合。(我们对量化文件进行了 gzip 压缩,使数据包更小,这对于我们用于导入文件的 R 函数来说不是问题。)

    samples<-read.table(file.path(dir,"samples.txt"),header = TRUE)
    files<-file.path(dir,"salmon",samples$run,"quant.sf.gz")
    names(files)<-paste0("samples",1:6)
    all(file.exists(files))
    

    在这里插入图片描述

    2、将转录本映射到基因

    转录本需要与基因 ID 相关联,以便进行基因级汇总。因此,我们将构建一个名为data.frame 的tx2gene两列:1) 转录本 ID 和 2) 基因 ID。列名无关紧要,但必须使用此列顺序。转录本 ID 必须与丰度文件中使用的相同。这可以通过在下载转录组 FASTA 的同时下载tx2geneGTF 文件并使用 Bioconductor 的TxDb基础设施从 GTF 文件生成来最容易地完成。

    创建tx2gene data.frame可以通过在TxDb对象上调用AnnotationDbi包中的select函数来完成。以下代码可用于构造这样的表:

    #将转录本映射到基因
    library("TxDb.Hsapiens.UCSC.hg38.knownGene")
    txdb<-TxDb.Hsapiens.UCSC.hg38.knownGene
    txdb
    k<-keys(txdb,keytype = "TXNAME")
    #k        创建tx2gene 
    tx2gene<-select(txdb,k,"GENEID","TXNAME")
    #使用Gencode V27 chr 转录本构建salmon索引
    library("readr")
    tx2gene<-read_csv(file.path(dir,"tx2gene.gencode.v27.csv"))
    head(tx2gene)
    

    3、tximport命令

    #tximport 命令——将salmon转录本量化导入R
    library("tximport")
    library("jsonlite")
    txi<-tximport(files,type="salmon",tx2gene = tx2gene)
    #txi 对象是一个矩阵列表
    names(txi)
    txi$counts[1:3,1:3]#提取前三行前三列的txi矩阵中的counts
    txi$length[1:3,1:3]#提取前三行前三列的txi矩阵中的length
    txi$abundance[1:3,1:3]#提取前三行前三列的txi矩阵中的abundance
    txi$countsFromAbundance
    

    在这里插入图片描述
    如果我们继续使用 GEUVADIS 样本,我们将使用以下代码行创建DESeqDataSet。因为样本之间没有差异(相同的群体和相同的测序批次),我们指定了一个~1的设计公式,这意味着我们只能拟合一个截距项——所以我们不能对这些样本进行差异表达分析。

    library("DESeq2")
    dds<-DESeqDataSetFromTximport(txi,samples,~1)
    dds$center
    dds$pop
    

    在这里插入图片描述
    告一段落,明天——>探索性数据分析

    展开全文
  • 生信入门(五)——使用DESeq2进行RNA-seq数据分析 文章目录生信入门(五)——使用DESeq2进行RNA-seq数据分析四、探索性数据分析五、差异数据分析六、AnnotationHub 本篇接上一篇,本篇做探索性数据分析,差异表达...

    生信入门(五)——使用DESeq2进行RNA-seq数据分析


    本篇接上一篇,本篇做探索性数据分析,差异表达分析以及后面步骤

    四、探索性数据分析

    1、简单EDA

    #探索性数据分析
    #BiocManager::install("airway")
    library("airway")
    data("airway")
    #指定untrt是dex 变量的参考级别
    airway$dex<-relevel(airway$dex,"untrt")
    airway$dex
    #快速检查与基因唯一对齐的数百万个片段(round的第二个参数告诉要保留多少小数点)
    round(colSums(assay(airway))/1e6,1)
    #通过拉出SummarizedExperiment的colData插槽来检查有关样本的信息:
    colData(airway)
    table(airway$cell)
    table(airway$dex)
    #加载DESeq2,创建DESeqDataSet,测试地塞米松治疗之间的差异
    library("DESeq2")
    dds<-DESeqDataSet(airway,design = ~cell+dex)
    #dds
    #执行最小过滤以减小数据集的大小
    #(如果4个或者更多样本的计数不为5或更多,我们不需要保留基因,因为这些基因将没有检测差异的统计能力,也没有计算样本之间距离的信息)
    keep<-rowSums(counts(dds)>=5)>=4
    table(keep)
    dds<-dds[keep,]
    #基本的探索性分析是检查每个样本计数的箱线图(对数据取对数,方便大基数主导箱线图)
    boxplot(log10(counts(dds)+1))
    dds<-estimateSizeFactors(dds)
    #归一化处理
    boxplot(log10(counts(dds,normalized=TRUE)+1))
    

    运行结果
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    2、EDA 的数据转换

    取计数的对数加上伪计数1是一种常见的变换,但是他往往会夸大低计数的采取方差,以至于它甚至大于样本组之间的生物变异。因此在DESeq2中,提供了产生对数尺度数据的转换,从而消除了系统趋势。一般推荐的转换是方差稳定转换或VST,他可以用VST函数调用:
    VST函数不返回DESeqDataSet,因为其不返回计数,而是返回连续值(在log2尺度上),可以使用分析访问转换后的数据。assay函数转换

    #EDA数据转换 (使用VST函数调用转换)
    vsd<-vst(dds)
    class(vsd)
    assay(vsd)[1:3,1:8]#使用VST函数返回的不是计数,而是返回连续值
    

    运行结果
    在这里插入图片描述

    3、主成分图

    VST数据适用于计算样本之间的距离或执行PCA(principal components analysis)。简而言之,PCA图使我们能够将数据中最主要的变异轴可视化,这对于质量控制和了解样本见差异在不同条件下和条件内的大小都有用。可以看到PC1(数据中变化的主轴)将处理过的和未处理过的样本分开:

    #生成主成分图
    plotPCA(vsd,"dex")
    

    运行结果
    在这里插入图片描述
    通过额外的ggplot2代码,可以指出哪些样本属于哪个细胞系。
    注意:不建议将转换后的数据用于主要差异表达分析。相反,我们将使用原始计数和广义线性模型(GLM),该模型考虑低计数和高计数的预期方差。

    # 增加额外GGPLOT2代码,指出哪些样本属于哪个细胞系
    library("ggplot2")
    pcaData<-plotPCA(vsd,intgroup=c("dex","cell"),returnData=TRUE)
    # pcaData
    percentVar<-round(100*attr(pcaData,"percentVar"))
    ggplot(pcaData,aes(x=PC1,y=PC2,color=dex,shape=cell))+geom_point(size=3)+
      xlab(paste0("PC1:",percentVar[1],"%variance"))+
      ylab(paste0("PC2:",percentVar[2],"%variance"))+
      coord_fixed()
    

    运行结果
    在这里插入图片描述

    五、差异数据分析

    1、标准DE步骤

    # 差异表达分析
    # 标准DE步骤
    dds<-DESeq(dds)
    res<-results(dds)
    # res包含每个基因的结果(与DESeqDataSet中的顺序相同)
    # 如果想查看顶级基因,可以这样排序
    head(res[order(res$pvalue),])
    # 绘制顶部基因的计数
    plotCounts(dds,which.min(res$pvalue),"dex")
    

    运行结果
    在这里插入图片描述
    在这里插入图片描述
    使用以下方法检查用于地塞米松治疗导致的所有log2倍数变化——LFC超过计数平均值plotMA

    # 使用以下方法检查用于地塞米松治疗导致的所有log2倍数变化(LFC)超过计数平均值plotMA
    plotMA(res,ylim=c(-5,5))
    

    运行结果
    在这里插入图片描述
    上面的MA绘图左侧有很多不显着的大型LFC(灰色点)。由于对数计数的不精确,这些点获得了很大的LFC 。为了获得更多信息可视化和精确的基因按效应大小排序(对数倍数变化有时称为效应大小)建议使用DESeq2的功能来缩小LFC。——方法apeglm收缩估计器,在DESeq2的lfcShrink函数中可用:

    # BiocManager::install("apeglm")
    library("apeglm")
    resultsNames(dds)
    res2<-lfcShrink(dds,coef="dex_trt_vs_untrt",type="apeglm")
    par(mfrow=c(1,2))
    plotMA(res,ylim=c(-3,3),main="No shrinkage")
    plotMA(res2,ylim=c(-3,3),main="apeglm")
    

    运行结果
    在这里插入图片描述

    2、最小效应量

    如果不想报告具有小LFC的重要基因,可以通过选择LFC并对此进行测试来指定最小的生物学意义效应大小。我们可以使用没有shrink的LFC或lfcShrink提供的LFC使用apeglm方法执行这样的阈值测试

    注意:测试针对LFC阈 ≠ 测试针对0假设然后过滤上LFC值
    apeglm方法提供S值时svalue=TRUE,或者当我们提供最小的影响大小(如上),这些类似于q值或者调整后的p值,因为s值 小于α应该有一个总的误号率或绝对值小于我们会给定的LFC阈值,其界限为α

    res.lfc<-results(dds,lfcThreshold = 1)
    res.lfc2<-lfcShrink(dds,coef = "dex_trt_vs_untrt",type="apeglm",lfcThreshold = 1)
    
    par(mfrow=c(1,2))
    plotMA(res.lfc,ylim=c(-5,5),main="NO shringkage ,LFC test")
    plotMA(res.lfc2,ylim=c(-5,5),main="apeglm,LFC test",alpha=0.01)
    

    在这里插入图片描述

    六、AnnotationHub

    1、查询AnnotationHub

    使用AnnotationHub包将附加信息附加到结果表中。AnnotationHub为超过40000条注释记录提供了一个易于使用的界面。一条记录可以来自ENCODE,人类基因组的序列,一个TXDB含有约转录物和基因,或信息OrgDb含有有关特定生物体的生物标识符的一般信息

    # annnotationhub
    # BiocManager::install("AnnotationHub",force = TRUE)
    library("AnnotationHub")
    ah<-AnnotationHub()
    

    使用以下代码启动浏览器,用于导航通过AnnotationHub可用的所有记录

    display(ah)
    

    运行结果
    在这里插入图片描述
    还可以使用带有查询功能的关键字进行查询

    query(ah,c("OrgDb","Homo sapiens"))
    

    运行结果
    在这里插入图片描述

    # 如果要下拉特定记录,可以使用双括号和记录名称
    hs<-ah[['AH92581']]
    hs
    

    运行结果
    在这里插入图片描述

    2、映射id

    # rowbames结果表是ENSEMEL IDs,大部分在OrgDb中(尽管有几千都不在)
    columns(hs)
    table(rownames(res)%in%keys(hs,"ENSEMBL"))
    #可以使用maplds函数添加基因符号,keytype=ENSEMBL,column="SYMBOL"
    res$symbol<-mapIds(hs,rownames(res),column = "SYMBOL",keytype = "ENSEMBL")
    head(res)
    

    运行结果
    在这里插入图片描述

    七、Building Report

    1、报告工具

    有许多从Bioconductor构建交互报告的软件包。其中两个是Reporting Tools和Glimma,他们都提供HTML报告,允许合作者检查基因组分析中的顶级基因(或任何感兴趣的特征)
    编译Reporting Tools报告的代码:

    # building report
    BiocManager::install("ReportingTools")
    library("ReportingTools")
    tmp<-tempdir()
    rep<-HTMLReport(shortName = "airway",title = "Airway DGE",
                    basePath = tmp,reportDirectory = "report")
    publish(res,rep,dds,n=20,make.plots=TRUE,factor=dds$dex)
    finish(rep)
    #在浏览器中启动
    browseURL(file.path(tmp,"report","airway.html"))
    

    运行结果

    在这里插入图片描述

    2、Glimma

    另一个可以生成交互式报告的包是Glimma。gIMDPlot 构建了一个交互式MA绘图,其中将鼠标悬停在左侧MA绘图中的基因上将显示右侧样本的计数。但即将在工具提示和屏幕底部的列表中显示基因信息。将鼠标悬停在右侧的样本上将在工具提示中提供样本ID

    # 还有Glimma生成交互式报告
    library("Glimma")
    status<-as.numeric(res$padj<.1)
    anno<-data.frame(GeneID=rownames(res),symbol=res$symbol)
    glMDPlot(res2,status=status,counts=counts(dds,normalized=TRUE),
             groups=dds$dex,transform=FALSE,
             samples=colnames(dds),anno=anno,
             path=tmp,folder="glimma",launch=TRUE)#当launch=FALSE时,不显示,可以结合下面语句显示,若为TRUE 则不需要下面语句
    #browseURL(file.path(tmp,"glimma","MD-Plot.html"))
    

    运行结果
    在这里插入图片描述

    八、与ZINB-WaVE的集成

    1、模拟飞溅

    在最后一小节,展示DESeq2与另一个Bioconductor包zinbwave集成,以便对额外的零进行建模和解释(超过负二项式模型的预期)——对于单细胞RNA-seq实验很有用
    此处,使用splatter包来模拟单细胞RNA-seq数据(Zappia,Phipson和Oshlack)
    **注:虽然 ZINB-WaVE 和 ZINGER 等方法可以成功识别多余的零点,但它们不能轻易区分其根本原因,即技术(例如,丢失)和生物(例如,爆裂)零点.
    **
    当感兴趣的信号不在零分量中时,可以使用下面概述的零通胀加权方法。也就是说,如果您想找到跨细胞群转录爆发的生物学差异,以下方法将无法帮助您找到这些差异。相反,它有助于发现除零分量之外的计数差异(这些零是生物的还是技术的)。

    2、用飞溅模拟单细胞计数数据

    # BiocManager::install("splatter")
    library("splatter")
    params<-newSplatParams()
    params<-setParam(params,"de.facLoc",1)
    params<-setParam(params,"de.facScale",.25)
    params<-setParam(params,"dropout.type","experiment")
    params<-setParam(params,"dropout.mid",3)
    set.seed(1)
    sim<-splatSimulate(params,group.prob=c(.5,.5),method = "groups")
    # sim
    plot(log10(rowMeans(assays(sim)[["TrueCounts"]])),
         rowMeans(assays(sim)[["Dropout"]]))
    

    运行结果
    在这里插入图片描述

    # 将存储真是的log2倍数变化进行比较
    rowData(sim)$log2FC<-with(rowData(sim),log2(DEFacGroup2/DEFacGroup1))
    
    # 负二项式分量在均值上的真实离差
    rowData(sim)$trueDisp<-rowMeans(assays(sim)[["BCV"]])^2
    gridlines<-c(1e-2,1e-1,1);cols<-c("blue","red","darkgreen")
    with(rowData(sim)[rowData(sim)$GeneMean>1,],
         plot(GeneMean,trueDisp,log="xy",xlim=c(1,300),ylim=c(.01,5)))
    abline(h=gridlines,col=cols)
    text(300,gridlines,labels=gridlines,col=cols,pos=3)
    

    运行结果
    在这里插入图片描述
    今日告一段落,回见。嘻嘻嘻。国庆放假了。祝大家假期快乐!

    展开全文
  • A survey of best practices for RNA-seq data analysis RNA-seq数据分析指南 ...A survey of best practices for RNA-seq data analysis,我把它叫做RNA-seq数据分析指南。这篇文章是由佛罗里达大学等单位的...
  • ##批处理比对数据文件 #!/bin/bash #$ -cwd #$ -j n #$ -S /bin/bash #$ -pe mpi 10 #$ -q NGS #$ -N gaomeiling #$ -e gaomeiling.err #$ -o gaomeiling.log wkpath=/gluster/home/gaomeiling/work/wancheng_liu/...
  • 生信入门(二)——使用limma、Glimma和edgeR,RNA-seq数据分析 文章目录生信入门(二)——使用limma、Glimma和edgeR,RNA-seq数据分析一、简介 一、简介 简单且高效地分析RNA 测序数据的能力是生信的核心优势之一。...
  • RNA-seq数据分析及其应用 克利夫兰诊所生物信息学研讨会#2(2018年7月) 这是我们在2018年举办的第二次生物信息学研讨会。这次我们将介绍基本的RNA序列分析。 没有任何先决条件。 对于已经为该研讨会注册的人员,...
  • RNA-Seq数据分析

    万次阅读 2014-10-13 15:09:32
    从原始的数据开始,进行reads回帖,到拼接转录本,计算表达量,分析差异表达,最后可视化分析结果。  TopHat是一个把reads回帖到基因组上的工具。首先用Bowtie把reads回帖到基因组上,然后通过拼接,我们就可以在...
  • 使用limma、Glimma和edgeR,RNA-seq数据分析 文章目录使用limma、Glimma和edgeR,RNA-seq数据分析六、差异表达分析1、创建设计矩阵和对比2、从表达计数数据中删除异方差3、拟合线性模型以进行比较4、检查DE基因数量5、...
  • 我开发了 scGEAToolbox——一个用于 scRNA-seq 数据分析的 Matlab 工具箱。 它包含用于数据归一化、特征选择、批量校正、插补、细胞聚类、轨迹/伪时间分析和网络构建的一整套功能,可以组合和集成来构建自定义工作...
  • 基因型组织表达项目的RNA-seq数据分析包含自闭症相关重复序列扩展的基因。 该项目的背景是一项遗传研究,由Trost B.等人在《自然》杂志2020年发表(Trost B,Engchuan W,Nguyen CM,Thiruvahindrapuram B,...
  • 使用limma、Glimma和edgeR对RNA-seq数据分析1. 摘要 Xueyi Dong^1, Charity Law^2, Monther Alhamdoosh3, Shian Su1, Luyi Tian2, Gordon K. Smyth4 and Matthew E. Ritchie5 1The Walter and Eliza Hall Institute ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,269
精华内容 907
关键字:

rna-seq数据分析