生信小白的RNA-seq实战历程

作者注：虽然是实战，其实只能称得上学习笔记而已，初学过程中参考了大量博客和帖子,还有大神引用的大牛的帖子,参考列表见最后。

1、软件安装

1.1 硬件系统情况

系统：BioLinux8（ubuntu14.04.2-x64）

吐槽个一下这个系统，普通帐户竟然环境变量错误（source ~/.bashrc)，对我等小白来说实在头大，只有用root用户运行了。

后来通过百度错误问题发现，这个系统用的是zsh，所以更新环境变量的代码是 source ~/.zshrc 这样就可以了。

电脑配置：

Lenovo-ThinkPad-E431

1. CPU系列 英特尔 酷睿i3 3代系列

2. CPU型号 Intel 酷睿i3 3120M

3. CPU主频 2.5GHz

4. 总线规格 DMI 5 GT/s

5. 三级缓存 3MB

6. 核心架构 Ivy Bridge

7. 核心/线程数 双核心/四线程

8. 制程工艺 22nm

9. 内存容量 12GB（4GB×1 + 8GB×1）

1.2 软件安装

有简单的不必去重新编译代码了，bioconda几个命令解决，可能软件依赖会有点问题。

1）Miniconda安装

1. wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh

2. bash Miniconda3-latest-Linux-x86_64.sh

3. source ~/.zshrc

一路Enter搞定

2）sratoolkit

conda install -c jfear sratoolkit

3）fastqc

conda install fastqc

发现fastqc是BioLinux自带的，执行这个命令只是升级了java

4）hisat2

conda install hisat2

hisat2 -h #测试

这有个小插曲，提示hisat2依赖于python3.5, 而我装的是3.6，于是百度得知用conda安装3.5版本的：

conda install python=3.5

5）samtools

1. conda install samtools

2. samtools

samtools也是自带的，执行这个命令会说依赖冲突，因为我的conda是Python3.6版，有点新，不兼容正常。

6）htseq-count

1. conda install -c bioconda htseq

2. #测试

3. python

4. >>> import HTSeq

7）R

BioLinux自带了，升级一下到3.4.1

1. sudo add-apt-repository ppa:marutter/rrutter

2. sudo apt-get update

3. sudo apt-get install r-base r-base-dev

8）rstudio

1. conda install rstudio

2. rstudio

---

2、读文章拿到测序数据

直接拿jimmy的脚本修改得来的。仔细分析项目数据后，得到要下载的数据是SRR3589958-62

这个脚本包含了下载数据，sra格式转化为fastq，fastqc对转化的数据进行分析，代码如下：

1. for ((i=958;i<=958;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP075/SRP075747/SRR3589$i/SRR3589$i.sra;done

2. ls *sra |while read id; do fastq-dump --gzip --split-3 $id;done

看到大神用的.gz格式才发现太省空间了，至少省了一半以上，膜拜！

---

3、了解fastq测序数据

1.fastqc

1. zcat SRR3589956_1.fastq.gz | fastqc -t 4 stdin

2. fastqc SRR3589956_1.fastq.gz

3. #t 线程，每个250M内存

2.multiQC

1. conda install multiqc

2. # 先获取QC结果

3. ls *gz | while read id; do fastqc -t 4 $id; done

4. # multiqc

5. multiqc *fastqc.zip --pdf

学习笔记是纸质版拍照的

---

4、了解参考基因组及基因注释

1）下载参考基因组

1. wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz #wget 下载或者其他工具下载

2. tar -zvf chromFa.tar.gz

3. cat *.fa > hg19.fa

4. rm chr*

2）下载基因组注释gtf文件

GTF和GFF之间的区别：

数据结构：都是由9列构成，分别是reference sequence name; annotation source; feature type; start coordinate; end coordinate; score; strand; frame; attributes.前8列都是相同的，第9列不同。

GFF第9列：都是以键值对的形式，键值之间用“=”连接，不同属性之间用“；”分隔，都是以ID这个属性开始。下图中有两个ID，说明是不同的序列。

GTF第9列：同样以键值对的形式，键值之间是以空格区分，值用双引号括起来；不同属性之间用“；”分隔；开头必须是geneid, transciptid两个属性。

1. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/GRCh37_mapping/gencode.v26lift37.annotation.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... 7.annotation.gtf.gz

3)下载IGV(Integrative Genomics Viewer)

下载地址为： http://software.broadinstitute.org/software/igv/download

4)genome -> Load Genome From Files加载之前得到基因组文件

5)Tool -> Run igvtools，进行排序

6)加载gff基因注释文件，File -> Load From Files

7)可视化分析 同样推荐阅读生信宝典公众号文章。

---

5、序列比对

5.1 下载index

1. cd referece && mkdir index && cd index

2. wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz

3. tar -zxvf hg19.tar.gz

5.2 比对

1. mkdir -p RNA-Seq/aligned

2. for i in seq ‘56 58’

3. do

4. hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam

5. done

5.3 HISAT2输出结果

5.4 格式转换，排序，索引

1. for i in `seq 56 58`

2. do

3.    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam

4.    samtools sort SRR35899${i}.bam SRR35899${i}_sorted.bam

5.    #这里我用的是0.1.19版本，不用加参数 -o    

6.    #ps:我加了-o之后重定向输出结果有5个G之工巨，xshell直接死机，直接运行，电脑终端一直不停跳乱码的东西。

7.    samtools index SRR35899${i}_sorted.bam

8. done

判断sam排序两种方式的不同：

1. head -100 SRR3589957.sam > test.sam

2. samtools view -b  test.sam > test.bam

3. samtools view test.bam | head

默认排序

1. samtools sort test.bam default

2. samtools view default.bam | head

Sort alignments by leftmost coordinates, or by read name when -n is used

1. samtools sort -n test.bam sort_left

2. samtools view sort_left.bam | head

5.5 samtools view

提取1号染色体1234-123456区域的比对read

samtools view SRR3589957_sorted.bam chr1:1234-123456 | head

flagstat看下总体情况

1. samtools flagstat SRR3589957_sorted.bam

用samtools筛选恰好配对的read,就需要用0x10

1. samtools view -b -f 0x10 SRR3589957_sorted.bam chr1:1234-123456  > flag.bam

2. samtools flagstat flag.bam

5.6 比对质控(QC)

1. # Python2.7环境下

2. pip install RSeQC

对bam文件进行统计分析(70~90的比对率要求)

1. bam_stat.py -i SRR3589956_sorted.bam

基因组覆盖率的QC

需要提供bed文件，可以直接RSeQC的网站下载(看文件只有1M多，就没有考虑转换):

1. wget https://downloads.sourceforge.net/project/rseqc/BED/Human_Ho<br>mo_sapiens/hg19_RefSeq.bed.gz

1. read_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed

5.7 IGV查看

载入参考序列，注释和BAM文件

---

6 reads计数

1. # 安装

2. conda install htseq

3. # 使用

4. # htseq-count [options] <alignment_file> <gtf_file>

5. htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count

循环处理多个BAM文件：

1. mkdir -p RNA-Seq/matrix/

2. for i in `seq 56 58`

3. do

4.    htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log

5. done

合并表达矩阵

在生信媛的python脚本上略做了修改

1. !/usr/bin/python3

2. def combine_count(count_list):

3.  mydict = {}

4.  for file in count_list:

5.      for line in open(file, 'r'):

6.          #print(line)

7.          if line.startswith('E'):

8.              key,value = line.strip().split('\t')

9.              if key in mydict:

10.                  mydict[key] = mydict[key] + '\t' + value

11.              else:

12.                  mydict[key] = value

13.  sorted_mydict = sorted(mydict)

14.  out = open('count_out.txt', 'w')

15.  for k in sorted_mydict:

16.      #print(k, mydict[k])

17.      #break

18.      out.write(k + '\t' + mydict[k] +'\n')

19. count_list = ['SRR3589956.count', 'SRR3589957.count', 'SRR3589958.count']

20. combine_count(count_list)

简单分析

这里由于对R接触还很少，不少命令不明白，只有运行一下试试

1) 导入数据

1. options(stringsAsFactors = FALSE)# import data if sample are small

2. control <- read.table("/media/biolinux/LENOVO_imcdul/biodata/SRR3589956.count",sep="\t", col.names = c("gene_id","control"))

2) 数据整合

1. # merge data and delete the unuseful row

2. raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")

3. raw_count_filt <- raw_count[-1:-5,]

4. ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)

5. row.names(raw_count_filt) <- ENSEMBL

3） 总体情况

1. summary(raw_count_filt)

4）看几个具体基因

1. > GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]

2. #EBI上搜索GAPDH找到ID为ENSG00000111640

3. > GAPDH

文章研究的AKAP95（ENSG00000105127）的表达量在KD中都降低了

1. > AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]

2. > AKAP95

---

参考内容列表

1、jimmy的导航贴 （http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost)

2、转录组（一）(http://www.biotrainee.com/thread-1800-1-1.html) ( HOPTOP )

3、转录组入门(1)(http://www.biotrainee.com/thread-1796-1-1.html)- (青山屋主)

4、转录组入门（1）Mac上软件准备作业 （http://www.biotrainee.com/thread-1810-1-1.html）

5、PANDA姐的转录组入门(1)：计算机资源的准备（http://www.biotrainee.com/thread-1847-1-1.html）

6、转录组作业（一）：来自零基础的小白（http://www.biotrainee.com/thread-1850-1-1.html）

7、转录组入门(2)：读文章拿到测序数据（http://www.biotrainee.com/thread-1743-1-1.html）

8、转录组入门(2)-（青山屋主)（http://www.biotrainee.com/thread-1884-1-1.html）

9、PANDA姐的转录组入门(2)：读文章拿到测序数据（http://www.biotrainee.com/thread-1853-1-1.html）

10、转录组入门(3)：了解fastq测序数据

11、转录组（三）：（HOPTOP）（http://www.biotrainee.com/thread-1831-1-1.html）

12、转录组入门(3)-（青山屋主）（http://www.biotrainee.com/thread-1885-1-1.html）

13、PANDA姐的转录组入门(3):了解fastq测序数据（http://www.biotrainee.com/thread-1853-1-1.html）

14、hoptop的：转录组作业（四）（http://www.biotrainee.com/thread-1835-1-1.html）

15、转录组入门(5)](http://www.biotrainee.com/thread-1870-1-1.html)： 序列比对（HOPTOP）

16、生信媛 转录组入门(6)： reads计数

17、http://www.bio-info-trainee.com/2218.html

---

猜你喜欢

菜鸟入门

数据分析

编程实践

直播基因组分析

生信技能树

编辑：思考问题的熊