RNAseq實(shí)際操作（實(shí)戰(zhàn)）
首先聲明击困，雖然是實(shí)戰(zhàn)，但是其實(shí)是學(xué)習(xí)筆記而已保礼，初學(xué)沛励，參考了大量大神的博客和帖子,還有大神引用的大牛的帖子,參考列表見最后。
1炮障、軟件安裝
1.1 硬件系統(tǒng)情況
系統(tǒng)：BioLinux8（ubuntu14.04.2-x64）
吐槽個(gè)一下這個(gè)系統(tǒng)目派，普通帳戶竟然環(huán)境變量錯(cuò)誤（source ~/.bashrc)，對(duì)我等小白來(lái)說(shuō)實(shí)在頭大胁赢，只有用root用戶運(yùn)行了企蹭。
后來(lái)通過(guò)百度錯(cuò)誤問題發(fā)現(xiàn)，這個(gè)系統(tǒng)用的是zsh，所以更新環(huán)境變量的代碼是 source ~/.zshrc 這樣就可以了谅摄。
電腦配置：
Lenovo-ThinkPad-E431
CPU系列 英特爾 酷睿i3 3代系列
CPU型號(hào) Intel 酷睿i3 3120M
CPU主頻 2.5GHz
總線規(guī)格 DMI 5 GT/s
三級(jí)緩存 3MB
核心架構(gòu) Ivy Bridge
核心/線程數(shù) 雙核心/四線程
制程工藝 22nm
內(nèi)存容量 12GB（4GB×1 + 8GB×1）
1.2 軟件安裝
有簡(jiǎn)單的不必去重新編譯代碼了徒河，bioconda幾個(gè)命令解決，可能軟件依賴會(huì)有點(diǎn)問題送漠。
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í)行這個(gè)命令只是升級(jí)了java
4）hisat2
conda install hisat2
hisat2 -h #測(cè)試
這有個(gè)小插曲，提示hisat2依賴于python3.5, 而我裝的是3.6闽寡，于是百度得知用conda安裝3.5版本的：
conda install python=3.5
5）samtools
conda install samtools
samtools
''#samtools也是自帶的代兵，執(zhí)行這個(gè)命令會(huì)說(shuō)依賴沖突，因?yàn)槲业腸onda是Python3.6版爷狈，有點(diǎn)新植影，不兼容正常。
6）htseq-count
conda install -c bioconda htseq
’#測(cè)試
python
>>> import HTSeq
7）R
BioLinux自帶了涎永，升級(jí)一下到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思币、讀文章拿到測(cè)序數(shù)據(jù)
直接拿jimmy的腳本修改得來(lái)的。仔細(xì)分析項(xiàng)目數(shù)據(jù)后羡微，得到要下載的數(shù)據(jù)是SRR3589958-62
這個(gè)腳本包含了下載數(shù)據(jù)谷饿，sra格式轉(zhuǎn)化為fastq，fastqc對(duì)轉(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測(cè)序數(shù)據(jù)（fastqc質(zhì)量控制學(xué)習(xí)）
1.fastqc
zcat SRR3589956_1.fastq.gz | fastqc -t 4 stdin
fastqc SRR3589956_1.fastq.gz
#t 線程锁荔，每個(gè)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列：都是以鍵值對(duì)的形式鸣峭，鍵值之間用“=”連接宏所，不同屬性之間用“；”分隔摊溶，都是以ID這個(gè)屬性開始爬骤。下圖中有兩個(gè)ID，說(shuō)明是不同的序列莫换。
GTF第9列：同樣以鍵值對(duì)的形式霞玄，鍵值之間是以空格區(qū)分骤铃，值用雙引號(hào)括起來(lái)；不同屬性之間用“坷剧；”分隔惰爬；開頭必須是<u>gene_id, transcipt_id</u>兩個(gè)屬性。
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惫企，進(jìn)行排序
6)加載gff基因注釋文件撕瞧，F(xiàn)ile -> Load From Files
7)可視化分析
同樣推薦閱讀生信寶典公眾號(hào)文章。
https://mp.weixin.qq.com/s/Q7pqycmQH58xU6hw_LECWA
5雅任、序列比對(duì)
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 比對(duì)
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é)果
兩個(gè)同時(shí)跑禽车，資源幾乎占滿。報(bào)錯(cuò)刊殉，以為硬件不足殉摔，后來(lái)才發(fā)現(xiàn)是硬盤爆了，于是上移動(dòng)硬盤记焊。
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è)G之工巨碗硬，xshell直接死機(jī)，直接運(yùn)行瓢颅，電腦終端一直不停跳亂碼的東西恩尾。
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號(hào)染色體1234-123456區(qū)域的比對(duì)read
samtools view SRR3589957_sorted.bam chr1:1234-123456 | head
flagstat看下總體情況
samtools flagstat SRR3589957_sorted.bam
用samtools篩選恰好配對(duì)的read,就需要用0x10
samtools view -b -f 0x10 SRR3589957_sorted.bam chr1:1234-123456 > flag.bam
samtools flagstat flag.bam
5.6 比對(duì)質(zhì)控(QC)
# Python2.7環(huán)境下
pip install RSeQC
對(duì)bam文件進(jìn)行統(tǒng)計(jì)分析(70~90的比對(duì)率要求)
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計(jì)數(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)處理多個(gè)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)
簡(jiǎn)單分析(這里由于對(duì)R接觸還很少，不少命令不明白信柿，只有運(yùn)行一下試試)
- 導(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"))
- 數(shù)據(jù)整合
# merge data and delete the unuseful row
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL
3） 總體情況
summary(raw_count_filt)
4）看幾個(gè)具體基因
> GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]
#EBI上搜索GAPDH找到ID為ENSG00000111640
> GAPDH
文章研究的AKAP95（ENSG00000105127）的表達(dá)量在KD中都降低了
> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> AKAP95
參考內(nèi)容列表：
1冀偶、http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost#jimmy的導(dǎo)航貼
2、轉(zhuǎn)錄組（一）](http://www.biotrainee.com/thread-1800-1-1.html) ( HOPTOP )
3渔嚷、轉(zhuǎn)錄組入門(1)](http://www.biotrainee.com/thread-1796-1-1.html)- (青山屋主)
4进鸠、轉(zhuǎn)錄組入門（1）Mac上軟件準(zhǔn)備作業(yè)
5、PANDA姐的轉(zhuǎn)錄組入門(1)：計(jì)算機(jī)資源的準(zhǔn)備
6圃伶、轉(zhuǎn)錄組作業(yè)（一）：來(lái)自零基礎(chǔ)的小白
7堤如、轉(zhuǎn)錄組入門(2)：讀文章拿到測(cè)序數(shù)據(jù)
[轉(zhuǎn)錄組入門（二） New(HOPTOP)
8蒲列、轉(zhuǎn)錄組入門(2)-（青山屋主)
9、PANDA姐的轉(zhuǎn)錄組入門(2)：讀文章拿到測(cè)序數(shù)據(jù)
10搀罢、轉(zhuǎn)錄組入門(3)：了解fastq測(cè)序數(shù)據(jù)
11蝗岖、轉(zhuǎn)錄組（三）：（HOPTOP）
12、轉(zhuǎn)錄組入門(3)-（青山屋主）
13榔至、PANDA姐的轉(zhuǎn)錄組入門(3):了解fastq測(cè)序數(shù)據(jù)
轉(zhuǎn)錄組入門(4)：了解參考基因組及基因注釋
14抵赢、hoptop的：轉(zhuǎn)錄組作業(yè)（四）
轉(zhuǎn)錄組入門(5)： 序列比對(duì)
15、轉(zhuǎn)錄組入門(5)](http://www.biotrainee.com/thread-1870-1-1.html)： 序列比對(duì)（HOPTOP）
轉(zhuǎn)錄組入門(6)： reads計(jì)數(shù)
16唧取、生信媛
轉(zhuǎn)錄組入門(7)： 差異基因分析
17铅鲤、生信媛
18、http://www.bio-info-trainee.com/2218.html
https://jiawen.zd200572.com
