计算生物学

计算生物学ex2

本实验使用Linux系统，华为云服务器，XShell连接
基因比对与文件操作
可视化使用Rstudio，IGV

EXERCISE1

PART 2

fastqc使用

cd FastQC
chmod 755 fastqc
export PATH=$PWD:$PATH

激活fastqc，可用which语句查看是否成功
用fastqc处理fastq文件，主要展示ERR45893一例

fastqc处理后生成了名为“ERR458493_fastqc.html”的文件
Xtpf7将文件下载至电脑进行查看

将文件下载到本地查看如图
一，显示基本统计信息

二，基本序列质量

三，序列得分
质量良好

其他情况：碱基含量正常；GC含量接近理论值，良好；其他指标正常，即质量良好。
##重复此步骤将余下fastq文件处理
后续实验将会用到，一开始未处理，后来回来重复操作，烦琐了很多

PART 3

ls -al

查看目录文件

解压查看文件头部语句如下

gunzip -c ERR458493.fastq.gz | head

生成结果如图所示

查看行数结果如下

root@ecs-kc1-large-2-linux-20201119153909:~/ex2# gunzip -c ERR458493.fastq.gz | wc -l
4375828

less查看R64文件

less R64.fa
less R64.gtf

空格符继续，q键退出查看
以>开头的为注释，后面是序列
fa文件展示如下

gtf文件展示
以#开头为注释

下一步：建立索引
mkdir genome #建立索引目录
STAR --runMode genomeGenerate #建索引
–runThreadN 32 #线程
–genomeDir genome #输出文件名
–genomeFastaFiles R64.fa #参考基因文件
–sjdbGTFfile R64.gtf #参考基因的注释文件
–sjdbOverhang 50 #值为read length-1
运行过程如下

下面进行比对

运行结果生成了新文件wt开头

less wt1_ReadsPerGene.out.tab语句查看文件，本例使用第二列

查看文件less wt1_Log.final.out

执行runSTAR.sh脚本将ERR文件全部处理

vim runSTAR.sh     #查看编辑脚本
：wq！     #保存并退出
bash runSTAR.sh  #执行脚本进行批处理

PART 4

用samtools生成索引文件

将文件下载到电脑
在IGV 创建基因数据库，操作如图
拷入bam文件（wt1，wt2，mu1）选择V区间，找寻三组基因的差异，本图片截取区间为V:14,131-22,783
由图可看出五号染色体YEL071W基因在wt表达较少，在mu表达较大

PART 5

脚本执行在PART3执行runSTAR.sh
合并count文件如下

root@ecs-kc1-large-2-linux-20201119153909:~/ex2# paste wt1_ReadsPerGene.out.tab wt2_ReadsPerGene.out.tab wt3_ReadsPerGene.out.tab mu1_ReadsPerGene.out.tab mu2_ReadsPerGene.out.tab mu3_ReadsPerGene.out.tab | \cut -f1,2,6,10,14,18,22 | \tail -n +5 > gene_count.txt

运行结果生成count文件

-rw-r--r--  1 root root    190783 Dec 14 21:13 gene_count.txt

将该文件下载至电脑打开生成的txt文件如图

PART6

由于服务器默认安装R环境为3.4不便于安装DEseq2包于是选择自己电脑Rstudio（4.0）进行实验
将gene_count.txt，和samples.txt文件下载至自己电脑，打开Rstudio运行
首先下载两个包，DEseq2安装方式可在csdn进行搜索。简单易懂，可用library（）语句查看安装是否成功
运行代码如下

library(DESeq2)
library(ggplot2)
cts <- as.matrix(read.csv("gene_count.txt",sep = "\t",row.names = 1,header = F))
coldata <- read.csv("samples.txt",sep = "\t",row.names = 1)
colnames(cts) <- rownames(coldata)
dds <- DESeqDataSetFromMatrix(countData = cts,colData = coldata,design = ~ Genotype)#转为D
vsd <- vst(dds,blind = F) #方差稳定变换
pcaData <- plotPCA(vsd,intgroup=c("Genotype"),returnData=T)
#intgroup: interesting groups: a character vector of names in colData(x
#returnData: should the function only return the data.frame of PC1 and
percentVar <- round(100*attr(pcaData,"percentVar"))
ggplot(pcaData,aes(PC1,PC2,color=Genotype))+
  geom_point(size=3)+
  xlim(-2.5,2.5)+ylim(-1,1)+
  xlab(paste0("PC1:",percentVar[1],"%variance"))+
  ylab(paste0("PC2:",percentVar[2],"%variance"))+
  geom_text(aes(label=name),vjust=2)

运行结果plot如图所示
之后建立基因表达差异表
代码如下

dds<-DESeq(dds)
resultsNames(dds) # lists the coefficients
# unfiltered results
res <- results(dds, name="Genotype_wt_vs_mu")
summary(res) # print the summary on screen
# filter the results for FDR < 0.05 and absolute-value-of-Fold-Change > 2
res05 <- res[which(res$padj<0.05 & abs(res$log2FoldChange)>log2(2)), ]
# sort by padj and write to file
res05Ordered <- res05[order(res05$padj),]
write.csv(as.data.frame(res05Ordered), file="wt_vs_mu_fdr05_fc2_results.csv")

查看导出csv文件，打开显示如图

FoldChange为差异倍数变化，即该值越大表示基因表达差异越明显，因此以此为依据降序排列

选择差异最为明显的前几基因id

在IGV可视化窗口中进行查看，在区间搜索直接搜索基因id即可
选择wt1和mu1进行比较
差异显著基因一SNF2（YOR290C）
该基因位于15号染色体区间XV:853,147-862,258

差异显著基因二PHO12（YHR215W）
该基因位于8号染色体区间VIII:550,691-554,910，同样是在wt表达量明显高于mu

表达基因可视化
包的安装同上
代码如下

res <- results(dds, name="Genotype_wt_vs_mu")
pdf("res_no_shrink.pdf")
plotMA(res, main = "No shrinkage", alpha=0.05, ylim=c(-4,4))
res_shrink <- lfcShrink(dds, coef="Genotype_wt_vs_mu", type="apeglm")
pdf("res_shrink.pdf")
plotMA(res_shrink, main = "Shrinkage by apeglm", alpha=0.05, ylim=c(-4,4))

生成pdf文档打开查看结果如下
不收缩的MA图
收缩的MA图

淘先锋技术网