生信分析-4

RNA-seq pipeline

数据分析流程

1. 原始数据
2. 数据质控
3. 比对质控数据
4. 比对上的数据的下游分析 – 基因表达差异分析、基因可变剪切、SNP/INDEL分析

软件安装(针对MAC系统)

用Anaconda快速搭建生物信息学分析平台

分析

• 下载fastq文件.fa

• 在终端输入以下命令

  

  raw _data 是fastq文件所在目录, 一般是压缩格式

  -q 是输出结果

  -t 4 是调用Mac电脑的4个核心跑数据, 注意自己电脑的CPU有几个核, 跑数据的核不能超过电脑的

  -o 表示输出结果到哪个目录, 后接目录路径

  ~/FastQC_result/raw/ 表示将结果输出到此目录下

  *.fq.gz 表示此目录下所有需要质控的fastq文件

• 输出的结果

  

• 查看质控报告 – 点击一个html文件

  

  注意: 这是其中一个图, 报告怎么看点击fastqc report

  1. 针对报告结果去除不正常的碱基

     $ fastx_trimmer [-h] [-f N] [-l N] [-z] [-v] [-i INFILE] [-o OUTFILE]
     [-h] = 获得帮助信息.
     [-f N] = 序列中从第几个碱基开始保留. 默认是1.
     [-l N] = 序列最后保留到多少个碱基，默认是整条序列全部保留.
     [-z] = 调用GZip软件，输出的文件自动经过压缩.

     zcat hela_ctrl_rep1_R1.fq.gz | fastx_trimmer -z -o hela_ctrl_rep1_R2_trimmer2.fq.gz -f 11 -l 140

     ==以上命令总是报错==

     该用下列方法

  • 先对.fa.gz文件解压缩

       gzip -d hela_ctrl_rep1_R2.fq.gz  # hela_ctrl_rep1_R2.fq.gz是文件名

  • 再用fasts_trimmer命令

      fastx_trimmer -i hela_ctrl_rep1_R2.fq -z  -o hela_ctrl_rep1_R2_trimmer.fq.gz -f 11 -l 140

  2. 去除reads两端的adapter

     cutadapt --times 1 -e 0.1 -0 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o hela_ctrl_rep1_R1_trimmer_adapt.fq.gz -p hela_ctrl_rep1_R2_trimmer_adapt.fq.gz hela_ctrl_rep1_R1_trimmer.fq.gz hela_ctrl_rep1_R2_trimmer.fq.gz

  报错

  

  *注意: * 上述-0 3命令写错了, 应该是大写的O

    cutadapt --times 1 -e 0.1 -O 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o hela_ctrl_rep1_R1_trimmer_adapt.fq.gz -p hela_ctrl_rep1_R2_trimmer_adapt.fq.gz hela_ctrl_rep1_R1_trimmer.fq.gz hela_ctrl_rep1_R2_trimmer.fq.gz

  –times 1 表示对每条reads进行一次cutadapt

  -e 0.1 表示允许adapt结果与测序序列相比有10%的错误率

  -O 3 表示adapt结果必须与测序序列有3个碱基以上的overlap

  –quality-cutoff 6 是常用指标

  -m 50 adapt后base 数低于50的丢弃

  -a AGATCGGAAGAGC 对应reads1 的通用引物

  -A AGATCGGAAGAGC 对应reads2 的通用引物

  -o 输出目录

• 比对rRNA index(就是rRNA的文件索引)

  • 比对不上的序列才能留下来 – 转录组

  • 关于如何下载参考基因组构成index文件, 请参考相关视频

  • 构建index命令

      # 先解压缩下载的fa.gz文件
      # 再合并解压缩后的文件
      cat chr1.fa chr2.fa ... > mm10_genome.fa
      # 建立index
      bowtie2-build ./mm10_genome.fa ./mm10_genome_bowtie2_index &  # ./mm10_genome.fa输入文件, ./mm10_genome_bowtie2_index输出文件

    • 比对

      nohup bowtie2 -x rRNA_index -1 hela_ctrl_rep1_R1_trimmer_adapt.fq.gz -2 hela_ctrl_rep1_R2_trimmer_adapt.fq.gz -S $sam_out -p 3 --un-conc-gz $fastq_unmap > $log 2>&1 &

      bowtie2 - x rRNA_index 是rRNA的索引文件, 也就是rRNA的参考序列

      -1 hela_ctrl_rep1_R1_trimmer_adapt.fq.gz 是reads1 trimmer后的结果

      -2 hela_ctrl_rep1_R2_trimmer_adapt.fq.gz 是reads2 trimmer后的结果

      -S $sam_out 表示比对结果的输出文件

      nohup表示在后台一直运行,即使退出程序

      -p 3 表示用3个cpu核心跑数据

      –un-conc-gz 表示比对不上rRNA的结果, 也就是转录组

      $fastq_unmap 将结果输出到此目录里

      *> $log 2>&1 & 将错误信息输出到日志文件里, 2表示标准错误流,错误信息, 1表示标准输出流,输出到屏幕 * 最后一个&表示后台运行

• 比对到参考基因组序列

    nohup tophat2 -p 3 -o output_dir/ $hg19_index un-conc-mate.1 un-conc-mate.2

  -p 3 表示用3个cpu核心跑数据

  -o output_dir 表示输出结果到此目录里

  $hg19_index 表示之前建立的基因组index文件

  un-conc-mate.1 un-conc-mate.2 表示没有比对上rRNA上的数据,也就是转录组数据

• 计算FPKM值

    cufflinks -o cufflinks_dir -p 3 -G /Users/zhangzhongnan/Downloads/RefSeq_genes_hg19.gtf ./*.bam

  -o cufflinks_dir 输出到此目录里

  -p 3 表示用3个cpu核心跑数据

  -G /Users/zhangzhongnan/Downloads/RefSeq_genes_hg19.gtf 表示转录组文件, 如何下载请Google

  ./*.bam表示之前比对后的结果

• 发现差异基因表达

  nohup cuffdiff -o $out_dir -p 4 --labels $label --min-reps-js-test 2 $hg19_gtf $ctrl_bam $treat_bam