【原】生信小白的RNA-seq實戰(zhàn)歷程

健明 2021-07-14

展開全文

作者注：雖然是實戰(zhàn)，其實只能稱得上學(xué)習(xí)筆記而已，初學(xué)過程中參考了大量博客和帖子,還有大神引用的大牛的帖子,參考列表見最后。

1、軟件安裝

1.1 硬件系統(tǒng)情況

系統(tǒng)：BioLinux8（ubuntu14.04.2-x64）

吐槽個一下這個系統(tǒng)，普通帳戶竟然環(huán)境變量錯誤（source ~/.bashrc)，對我等小白來說實在頭大，只有用root用戶運行了。

后來通過百度錯誤問題發(fā)現(xiàn)，這個系統(tǒng)用的是zsh，所以更新環(huán)境變量的代碼是 source ~/.zshrc 這樣就可以了。

電腦配置：

Lenovo-ThinkPad-E431

CPU系列英特爾酷睿i3 3代系列
CPU型號 Intel 酷睿i3 3120M
CPU主頻 2.5GHz
總線規(guī)格 DMI 5 GT/s
三級緩存 3MB
核心架構(gòu) Ivy Bridge
核心/線程數(shù) 雙核心/四線程
制程工藝 22nm
內(nèi)存容量 12GB（4GB×1 + 8GB×1）

1.2 軟件安裝

有簡單的不必去重新編譯代碼了，bioconda幾個命令解決，可能軟件依賴會有點問題。

1）Miniconda安裝

wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.zshrc

一路Enter搞定

2）sratoolkit

conda install -c jfear sratoolkit

3）fastqc

conda install fastqc

發(fā)現(xiàn)fastqc是BioLinux自帶的，執(zhí)行這個命令只是升級了java

4）hisat2

conda install hisat2

hisat2 -h #測試

這有個小插曲，提示hisat2依賴于python3.5, 而我裝的是3.6，于是百度得知用conda安裝3.5版本的：

conda install python=3.5

5）samtools

conda install samtools
samtools

samtools也是自帶的，執(zhí)行這個命令會說依賴沖突，因為我的conda是Python3.6版，有點新，不兼容正常。

6）htseq-count

conda install -c bioconda htseq
#測試
python
>>> import HTSeq

7）R

BioLinux自帶了，升級一下到3.4.1

sudo add-apt-repository ppa:marutter/rrutter
sudo apt-get update
sudo apt-get install r-base r-base-dev

8）rstudio

conda install rstudio
rstudio

2、讀文章拿到測序數(shù)據(jù)

直接拿jimmy的腳本修改得來的。仔細(xì)分析項目數(shù)據(jù)后，得到要下載的數(shù)據(jù)是SRR3589958-62

這個腳本包含了下載數(shù)據(jù)，sra格式轉(zhuǎn)化為fastq，fastqc對轉(zhuǎn)化的數(shù)據(jù)進(jìn)行分析，代碼如下：

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
ls *sra |while read id; do fastq-dump --gzip --split-3 $id;done

看到大神用的.gz格式才發(fā)現(xiàn)太省空間了，至少省了一半以上，膜拜！

3、了解fastq測序數(shù)據(jù)

1.fastqc

zcat SRR3589956_1.fastq.gz | fastqc -t 4 stdin
fastqc SRR3589956_1.fastq.gz
#t 線程，每個250M內(nèi)存

2.multiQC

conda install multiqc
# 先獲取QC結(jié)果
ls *gz | while read id; do fastqc -t 4 $id; done
# multiqc
multiqc *fastqc.zip --pdf

學(xué)習(xí)筆記是紙質(zhì)版拍照的

4、了解參考基因組及基因注釋

1）下載參考基因組

wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz #wget 下載或者其他工具下載
tar -zvf chromFa.tar.gz
cat *.fa > hg19.fa
rm chr*

2）下載基因組注釋gtf文件

GTF和GFF之間的區(qū)別：
數(shù)據(jù)結(jié)構(gòu)：都是由9列構(gòu)成，分別是reference sequence name; annotation source; feature type; start coordinate; end coordinate; score; strand; frame; attributes.前8列都是相同的，第9列不同。
GFF第9列：都是以鍵值對的形式，鍵值之間用“=”連接，不同屬性之間用“；”分隔，都是以ID這個屬性開始。下圖中有兩個ID，說明是不同的序列。
GTF第9列：同樣以鍵值對的形式，鍵值之間是以空格區(qū)分，值用雙引號括起來；不同屬性之間用“；”分隔；開頭必須是geneid, transciptid兩個屬性。

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./software/igv/download

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

5)Tool -> Run igvtools，進(jìn)行排序

6)加載gff基因注釋文件，F(xiàn)ile -> Load From Files

7)可視化分析同樣推薦閱讀生信寶典公眾號文章。

https://mp.weixin.qq.com/s/Q7pqycmQH58xU6hw_LECWA

5、序列比對

5.1 下載index

cd referece && mkdir index && cd index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar -zxvf hg19.tar.gz

5.2 比對

mkdir -p RNA-Seq/aligned
for i in seq '56 58’
do
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
done

5.3 HISAT2輸出結(jié)果

5.4 格式轉(zhuǎn)換，排序，索引

for i in `seq 56 58`
do
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
samtools sort SRR35899${i}.bam SRR35899${i}_sorted.bam
#這里我用的是0.1.19版本，不用加參數(shù) -o
#ps:我加了-o之后重定向輸出結(jié)果有5個G之工巨，xshell直接死機(jī)，直接運行，電腦終端一直不停跳亂碼的東西。
samtools index SRR35899${i}_sorted.bam
done

判斷sam排序兩種方式的不同：

head -100 SRR3589957.sam > test.sam
samtools view -b test.sam > test.bam
samtools view test.bam | head

默認(rèn)排序

samtools sort test.bam default
samtools view default.bam | head

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

samtools sort -n test.bam sort_left
samtools view sort_left.bam | head

5.5 samtools view

提取1號染色體1234-123456區(qū)域的比對read

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

flagstat看下總體情況

samtools flagstat SRR3589957_sorted.bam

用samtools篩選恰好配對的read,就需要用0x10

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

5.6 比對質(zhì)控(QC)

# Python2.7環(huán)境下
pip install RSeQC

對bam文件進(jìn)行統(tǒng)計分析(70~90的比對率要求)

bam_stat.py -i SRR3589956_sorted.bam

基因組覆蓋率的QC

需要提供bed文件，可以直接RSeQC的網(wǎng)站下載(看文件只有1M多，就沒有考慮轉(zhuǎn)換):

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

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

5.7 IGV查看

載入?yún)⒖夹蛄校⑨尯虰AM文件

6 reads計數(shù)

# 安裝
conda install htseq
# 使用
# htseq-count [options] <alignment_file> <gtf_file>
htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count

循環(huán)處理多個BAM文件：

mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
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
done

合并表達(dá)矩陣

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

!/usr/bin/python3
def combine_count(count_list):
mydict = {}
for file in count_list:
for line in open(file, 'r'):
#print(line)
if line.startswith('E'):
key,value = line.strip().split('\t')
if key in mydict:
mydict[key] = mydict[key] + '\t' + value
else:
mydict[key] = value
sorted_mydict = sorted(mydict)
out = open('count_out.txt', 'w')
for k in sorted_mydict:
#print(k, mydict[k])
#break
out.write(k + '\t' + mydict[k] +'\n')
count_list = ['SRR3589956.count', 'SRR3589957.count', 'SRR3589958.count']
combine_count(count_list)

簡單分析

這里由于對R接觸還很少，不少命令不明白，只有運行一下試試

1) 導(dǎo)入數(shù)據(jù)

options(stringsAsFactors = FALSE)# import data if sample are small
control <- read.table("/media/biolinux/LENOVO_imcdul/biodata/SRR3589956.count",sep="\t", col.names = c("gene_id","control"))

2) 數(shù)據(jù)整合