1. 斑馬魚基因組文件文件下載

ftp://ftp.ensembl.org/pub/release-97/fasta/danio_rerio/dna/

image.png

image.png

斑馬魚注釋文件
[ftp://ftp.ensembl.org/pub/release-97/gtf/danio_rerio/)

image.png

參考基因組下載 有三大全文網(wǎng)站提供參考基因組下載，它們分別是：
1.NCBI （https://www.ncbi.nlm.nih.gov/grc）
2.UCSC (http://hgdownload.soe.ucsc.edu/downloads.html)
3.Ensemble （http://asia.ensembl.org/index.html?redirect=no）
ps:目前最常用的人和小鼠的參考基因組版本如下：
|NCBI | UCSC| Ensemble| |GRCh36 | hg18 | ENSEMBL release_52 | |GRCh37 | hg19 | ENSEMBL release_59/61/64/68/69/75| |GRCh38 | hg38 | ENSEMBL release_76/77/78/80/81/82|

gz 格式文件，解壓縮命令是 $ gunzip 文件名

- 安裝bowtie2

$ sudo -s
##輸入密碼
$ sudo wget https://jaist.dl.sourceforge.net/project/bowtie-bio/bowtie2/2.3.4.1/bowtie2-2.3.4.1-linux-x86_64.zip
$ unzip bowtie2-2.3.4.1-linux-x86_64.zip

$ cd ~ #該環(huán)境變量，先回到主頁
$ vi .bashrc   

也可以不返回主頁進(jìn)入vi （$ vi ~/.bashrc） 末尾加上
-- 大寫G進(jìn)入最末端，小寫i進(jìn)入編輯模式，或者按兩次i。
按ESC，保存退出:wq ，不保存退出:q！

圖片4.png

如果出現(xiàn)上圖這種錯(cuò)誤，那是因?yàn)槎嗔藄wap文件
$ rm .bashrc.swp #刪除即可

圖片5.png

export PATH="$PATH:/home/qinghai/biosoft/bowtie2/bowtie2-2.3.5.1-linux-x86_64/"

操作 （一）進(jìn)入 1.vim filename （二）退出 1.:wq 末行模式，保存退出 2.:q 末行模式，直接退出 3.:q! 末行模式，不保存，強(qiáng)制退出 （三）輸入模式（在命令模式下操作） 1. i 從光標(biāo)所在位置前面開始插入 2.I 在當(dāng)前行首插入 3.a 從光標(biāo)所在位置后面開始輸入 4.A 在當(dāng)前行尾插入 5.o 在光標(biāo)所在行下方新增一行并進(jìn)入輸入模式 6.O 在當(dāng)前上面一行插入 （四）移動(dòng)光標(biāo)（在命令模式下操作） 1. gg 到文件第一行 2. G 到文件最后一行 (Shift + g) 3.^ 非空格行首 4.0 行首(數(shù)字0) 5.PATH:/home/qinghai/zebra/bowtie2
按Esc 退出編輯模式，輸入:wq 退出并保存

$ source ~/.bashrc
$ bowtie2 ##檢測(cè)是否安裝成功

• 輸入 echo $PATH來查看你的環(huán)境變量
  可以查看當(dāng)前的環(huán)境變量

 echo $PATH

export PATH=“$PATH:/home/qinghai/packages/bowtie2/bowtie2-2.3.4.1-linux-x86_64/”

$ mkdir biosoft  ##當(dāng)前路徑下建立文件夾
$ cd biosoft  ##進(jìn)入創(chuàng)建的文件夾下

另外一種方法安裝bowtie2 和samtools

采用系統(tǒng)自帶的 sudo apt install bowtie2
sudo apt install samtools

2. 建立索引

根據(jù)上面提示下載基因組序列（MT）

$ bowtie2-build Danio_rerio.GRCz11.dna.chromosome.MT.fa.gz genome
或者
$ bowtie2-build ./zebra_gtf/Danio_rerio_Ensemble_97.fa.gz genome

• 單獨(dú)線粒體序列生成索引比較快，全基因組比較慢
  生成包含6個(gè)文件的索引
  1/2/3/4 bt2 rev.1/2.bt2

3. bowtie2 進(jìn)行比對(duì)，得到sam文件

bowtie2 -p 6 -3 5 --local -x genome -1 MCK_1.clean.fq -2 MCK_2.clean.fq -S MCK.sam 

測(cè)序公司返回來數(shù)據(jù)，讀取結(jié)果如下

圖片3.png

• 測(cè)序公司返回的雙端測(cè)序原始clean Data（MCK_1.clean.fq MCK_2.clean.fq ）得到一個(gè)MCK.sam 文件，
• 多個(gè)樣的測(cè)序數(shù)據(jù)，可以批量文件進(jìn)行處理。也可以使用管道命令|直接將跳過生成dam文件，直接原型2個(gè)程序，直接得到較小文件bam。操作見下面

for ele in {44..55}#### for 循環(huán)，{}代表變化的集合
do 
bowtie2 -p 6 -3 5 --local -x genome  -1 SRR15524$ele.fastq.gz -2 SRR15524$ele.fastq.gz  –S example$ele.sam  #$變量替代相應(yīng)的數(shù)字
done# 運(yùn)行循環(huán)
for ele in {44..55}    ####批量將sam文件變?yōu)锽AM文件
do 
samtools sort example$ele.sam>./bam/example$ele.bam  ###批量將文件放進(jìn)一個(gè)文件夾，需要提前建立文件夾。
done

ls *.bam | while read id; do samtools index $id; done

############for ele in {MT1.MT2,MT3,WT1,WT2,WT3} do bowtie2 -p 6 -3 5 --local -x genome -1 ./data/M_H_ele_2.fq.gz | samtools sort O bam -@ 10 -o - > MTO1$ele.bam done

bowtie2 知識(shí) bowtie2與samtools兩個(gè)同時(shí)運(yùn)行，生成bam文件
單末端測(cè)序
$ bowtie2 -p 10 -x genome_index -U input.fq | samtools sort -O bam -@ 10 -o - > output.bam

雙末端測(cè)序
$ bowtie2 -p 10 -x genome_index -1 input_1.fq -2 input_2.fq | samtools sort -O bam -@ 10 -o - > output.bam

‘這條命令把bowtie2 生成的sam文件通過管道|傳遞到samtools，再將sam轉(zhuǎn)換為bam文件，省去中間sam文件的空間占genome_index 指的是用于bowtie2的索引文件。而不是參考基因組本身，構(gòu)建過程參考后文。genome_index 需要指定路徑及其共用文件名，比如我的索引文件放在/data/ref/bowtie2/mm10目錄下，但是需要輸入的參數(shù)

samtools 安裝

 $ wget https://sourceforge.net/projects/samtools/files/samtools/1.5/samtools-1.5.tar.bz2/download
$ tar -jxf download

export PATH="PATH:/home/qinghai/biosoft/samtools-1.5/"

samtools 將sam文件變?yōu)閎am文件

$ samtools sort MCK.sam > example.bam

IGV可視化匹配好的bam文件

另外還需要構(gòu)建索引，并和bam文件放在一塊才可以被IGV軟件打開

$ samtools index example.bam 
或samtools index MCK.bam

得到兩個(gè)文件，bam文件和bam.bai 文件

圖片2.png

這步不要4. count ---stringtie 計(jì)數(shù)，直接到featureCounts

安裝stringtie 軟件

mkdir stringtie-1.3-source   在已有的目錄下創(chuàng)建文件夾
wget http://ccb.jhu.edu/software/stringtie/dl/stringtie-1.3.3b.Linux_x86_64.tar.gz
tar zxvf stringtie-1.3.3b.Linux_x86_64.tar.gz

下載安裝包，解壓，記住安裝路徑，放在bashrc文件的最后面, 參考上面的操作
export PATH=$PATH:/home/qinghai/zebra/stringtie

• 進(jìn)行stringtie分析

 $ stringtie example.bam -p 16 -G Danio_rerio_Ensemble_97.gtf -B -o ./A/TRANSCK.gtf
或
 $ stringtie exampleHOM.bam -p 16 -G Danio_rerio_Ensemble_97.gtf -B -o ./A/TRANHOM.gtf
或
$ stringtie MCK.bam -p 16 -G ./GTF/Danio_rerio_Ensemble_97.gtf.gz -B -o ./a/MCK.gtf  
 ##所以文件需要先解壓，要不然會(huì)出錯(cuò)

圖片3.png

這里生成了結(jié)果gtf文件和ballgown需要的.ctab文件，還有基因的表達(dá)量文件gene_abund.tab，該文件包括基因的表達(dá)量FPKM以及TPM等。當(dāng)然如果你想要轉(zhuǎn)錄本的表達(dá)量，直接打開t_data.ctab這個(gè)文件，這里面有轉(zhuǎn)錄本的FPKM值。

將2個(gè)gtf文件進(jìn)行合并(或多個(gè)), 差不多如下面的格式

圖片1.png

$ stringtie --merge -p 8 -G Danio_rerio_Ensemble_97.gtf -o ./A/stringtie_merged.gtf ./A/mergelist.txt.txt
或
stringtie --merge -p 8 -G ./GTF/Danio_rerio_Ensemble_97.gtf -o ./A/stringtie_merged_MCK_HOM.gtf ./a/mergelist.txt

- 參考之前的內(nèi)容批量做

stringtie ./bam/example45.bam -p 16 -G ./A/MUS.gtf -B -o TRANS.gtf

for ele in {MT1,MT2,MT3,WT1,WT2,WT3}; do stringtie -p 10 -G Danio_rerio_Ensemble_97.gtf -B -o MH_$ele.gtf -l MH_$ele ./bam/M_H_$ele.bam; done 

stringtie -p 10 -G /ref/Mus_musculus/UCSC/mm10/Annotation/Genes/genes.gtf -o [output_gtf_file] -l [filename] [input_bam_file]
stringtie -p 8 -G ./genes/chrX.gtf -o ERR${i}_chrX.gtf -l ERR${i} ERR${i}_chrX.bam

用ballow或

$ for i in {HOM,SCK}; do stringtie -e -B -p -8 -G ./A/stringtie_merged.gtf -o ballgown/TRAN${i}/TRAB${i}_chrX.gtf example${i}.bam; done

4. featureCounts。輸出read是，txt格式，每個(gè)文件進(jìn)行counts

for ele in {MCK,MHOM}; do featureCounts -p -T 6 -t exon -g gene_id -a ./GTF/Danio_rerio_Ensemble_97.gtf -o ./count/counts$ele.txt $ele.bam; done 
或
for ele in {HOM,SCK}; do featureCounts -p -T 6 -t exon -g gene_id -a Danio_rerio_Ensemble_97.gtf -o ./B/counts$ele.txt example$ele.bam; done   #####每個(gè)文件進(jìn)行counts

所有輸入文件匯成一個(gè)txt

featureCounts -p -T 6 -t exon -g gene_id -a ./GTF/Danio_rerio_Ensemble_97.gtf -o ./count/counts_all.txt  *.bam 

批量進(jìn)行

for ele in {HOM,SCK}; do featureCounts -p -T 6 -t exon -g gene_id -a Danio_rerio_Ensemble_97.gtf -o ./B/counts_all.txt example$ele.bam; done   ###所有輸入文件匯成一個(gè)txt

最終結(jié)果

圖片4.png

4. 轉(zhuǎn)到R語言上

安裝DESeq2包 start R (version "3.6") and enter

>if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages

>("BiocManager")  BiocManager::install("DESeq2")

改變當(dāng)前路徑

> getwd()    #獲取當(dāng)前路徑
[1] "C:/Users/zqh19/Documents"
> setwd("E:/R/example")  #改變路徑
> getwd()
[1] "E:/R/example

圖片5.png

讀取txt 或 csv文件

data1<-read.table("counts_all.txt")
或
data1<-read.csv("counts_all.csv")
head(data1)

----dim(data1)###查看行數(shù)列數(shù) >head(data1) ##前六行數(shù)據(jù)

mycounts<-data1[,-(2:6)] ##先去掉第2-6列
write.csv(mycounts, "results2.csv",row.names=FALSE)
#重新寫入csv格式中，可手動(dòng)打開，去除行名，并將行名改為樣品名等，并刪除第一行數(shù)據(jù)
rownames(mycounts)<-mycounts[,1]   ##第一列數(shù)據(jù)命名為行名，并刪除第一列
mycounts<-mycounts[,-1]

head(mycounts) 或write.csv(mycounts, "results3.csv",row.names=FALSE) 查看數(shù)據(jù)的結(jié)果

R 軟件包下載

1.安裝Bioconductor

>if (!requireNamespace("BiocManager", quietly = TRUE))
>options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")  
>install.packages("BiocManager")
>library("BiocManager")

2.使用Bioconductor安裝org.Hs.eg.db：

>BiocManager::install("org.Hs.eg.db")

> getwd()
[1] "E:/R/example"
> data1<-read.table("counts_all.txt")
> mycounts<-data1[,-(2:6)]
> rownames(mycounts)<-mycounts[,1]
> inputdata<-as.data.frame(row.names(mycounts))
> head(inputdata)
  row.names(mycounts)
1              Geneid
2  ENSDARG00000099104
3  ENSDARG00000102407
4  ENSDARG00000102097
5  ENSDARG00000099319
6  ENSDARG00000099640
> newinput <- as.data.frame(gsub("\\.\\d*","",inputdata[,1]))
> my_dataset<-useDataset("drerio_gene_ensemble",mart=my_mart)

The downloaded source packages are in
‘C:\Users\zqh19\AppData\Local\Temp\RtmpgPdPdZ\downloaded_packages’

5. GENEid轉(zhuǎn)換 R

library('biomaRt')
library("curl")
my_mart <-useMart("ensembl")
datasets <- listDatasets(my_mart)
View(datasets)##查看數(shù)據(jù)庫

ensemble數(shù)據(jù)庫中包含有77個(gè)數(shù)據(jù)庫

圖片7.png

斑馬魚數(shù)據(jù)庫如下

圖片8.png

my_dataset <-useDataset("drerio_gene_ensembl",mart=my_mart)

根據(jù)基因ID獲取基因名，描述信息等

my_newid<- getBM(attributes = c("ensembl_gene_id","external_gene_name","description"),filters = "ensembl_gene_id",values=newinput,mart=my_dataset)

圖片9.png

查看數(shù)據(jù)庫的Attributes

listAttributes(my_dataset)

圖片10.png

兩個(gè)數(shù)據(jù)的merge （my_newid，，mycounts）

圖片1.png

final_res <- merge(my_newid,mycounts,by.x='ensembl_gene_id',by.y='V1')
head(final_res)
或
 write.csv(final_res,'A.csv')##寫入csv中查看

圖片2.png

注意control要放在前面####

inputdata<-as.data.frame(row.names(mycounts))
condition <- factor(c(rep("MCK",1),rep("MHOM",1)), levels = c("MCK","MHOM"))
my_dataset<-useDataset("drerio_gene_ensemble",mart=my_mart)
condition <- factor(c(rep("MCK",1),rep("MHOM",1)), levels = c("MCK","MHOM"))

dds <- DESeqDataSetFromMatrix(mycounts, inputdata, design= ~ condition)