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

参考基因组下载 有三大全文网站提供参考基因组下载，它们分别是：
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 ~ #该环境变量，先回到主页
$ vi .bashrc   

也可以不返回主页进入vi （$ vi ~/.bashrc） 末尾加上
-- 大写G进入最末端，小写i进入编辑模式，或者按两次i。
按ESC，保存退出:wq ，不保存退出:q！

图片4.png

如果出现上图这种错误，那是因为多了swap文件
$ rm .bashrc.swp #删除即可

图片5.png

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

操作 （一）进入 1.vim filename （二）退出 1.:wq 末行模式，保存退出 2.:q 末行模式，直接退出 3.:q! 末行模式，不保存，强制退出 （三）输入模式（在命令模式下操作） 1. i 从光标所在位置前面开始插入 2.I 在当前行首插入 3.a 从光标所在位置后面开始输入 4.A 在当前行尾插入 5.o 在光标所在行下方新增一行并进入输入模式 6.O 在当前上面一行插入 （四）移动光标（在命令模式下操作） 1. gg 到文件第一行 2. G 到文件最后一行 (Shift + g) 3.^ 非空格行首 4.0 行首(数字0) 5.PATH:/home/qinghai/zebra/bowtie2
按Esc 退出编辑模式，输入:wq 退出并保存

$ source ~/.bashrc
$ bowtie2 ##检测是否安装成功

• 输入 echo $PATH来查看你的环境变量
  可以查看当前的环境变量

 echo $PATH

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

$ mkdir biosoft  ##当前路径下建立文件夹
$ cd biosoft  ##进入创建的文件夹下

另外一种方法安装bowtie2 和samtools

采用系统自带的 sudo apt install bowtie2
sudo apt install samtools

2. 建立索引

根据上面提示下载基因组序列（MT）

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

• 单独线粒体序列生成索引比较快，全基因组比较慢
  生成包含6个文件的索引
  1/2/3/4 bt2 rev.1/2.bt2

3. bowtie2 进行比对，得到sam文件

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

测序公司返回来数据，读取结果如下

图片3.png

• 测序公司返回的双端测序原始clean Data（MCK_1.clean.fq MCK_2.clean.fq ）得到一个MCK.sam 文件，
• 多个样的测序数据，可以批量文件进行处理。也可以使用管道命令|直接将跳过生成dam文件，直接原型2个程序，直接得到较小文件bam。操作见下面

for ele in {44..55}#### for 循环，{}代表变化的集合
do 
bowtie2 -p 6 -3 5 --local -x genome  -1 SRR15524$ele.fastq.gz -2 SRR15524$ele.fastq.gz  –S example$ele.sam  #$变量替代相应的数字
done# 运行循环
for ele in {44..55}    ####批量将sam文件变为BAM文件
do 
samtools sort example$ele.sam>./bam/example$ele.bam  ###批量将文件放进一个文件夹，需要提前建立文件夹。
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 知识 bowtie2与samtools两个同时运行，生成bam文件
单末端测序
$ bowtie2 -p 10 -x genome_index -U input.fq | samtools sort -O bam -@ 10 -o - > output.bam

双末端测序
$ 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转换为bam文件，省去中间sam文件的空间占genome_index 指的是用于bowtie2的索引文件。而不是参考基因组本身，构建过程参考后文。genome_index 需要指定路径及其共用文件名，比如我的索引文件放在/data/ref/bowtie2/mm10目录下，但是需要输入的参数

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文件变为bam文件

$ samtools sort MCK.sam > example.bam

IGV可视化匹配好的bam文件

另外还需要构建索引，并和bam文件放在一块才可以被IGV软件打开

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

得到两个文件，bam文件和bam.bai 文件

图片2.png

这步不要4. count ---stringtie 计数，直接到featureCounts

安装stringtie 软件

mkdir stringtie-1.3-source   在已有的目录下创建文件夹
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

• 进行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  
 ##所以文件需要先解压，要不然会出错

图片3.png

这里生成了结果gtf文件和ballgown需要的.ctab文件，还有基因的表达量文件gene_abund.tab，该文件包括基因的表达量FPKM以及TPM等。当然如果你想要转录本的表达量，直接打开t_data.ctab这个文件，这里面有转录本的FPKM值。

将2个gtf文件进行合并(或多个), 差不多如下面的格式

图片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

- 参考之前的内容批量做

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格式，每个文件进行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   #####每个文件进行counts

所有输入文件汇成一个txt

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

批量进行

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   ###所有输入文件汇成一个txt

最终结果

图片4.png

4. 转到R语言上

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

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

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

改变当前路径

> getwd()    #获取当前路径
[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)###查看行数列数 >head(data1) ##前六行数据

mycounts<-data1[,-(2:6)] ##先去掉第2-6列
write.csv(mycounts, "results2.csv",row.names=FALSE)
#重新写入csv格式中，可手动打开，去除行名，并将行名改为样品名等，并删除第一行数据
rownames(mycounts)<-mycounts[,1]   ##第一列数据命名为行名，并删除第一列
mycounts<-mycounts[,-1]

head(mycounts) 或write.csv(mycounts, "results3.csv",row.names=FALSE) 查看数据的结果

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转换 R

library('biomaRt')
library("curl")
my_mart <-useMart("ensembl")
datasets <- listDatasets(my_mart)
View(datasets)##查看数据库

ensemble数据库中包含有77个数据库

图片7.png

斑马鱼数据库如下

图片8.png

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

根据基因ID获取基因名，描述信息等

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

图片9.png

查看数据库的Attributes

listAttributes(my_dataset)

图片10.png

两个数据的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)