其他分享
首页 > 其他分享> > 计算生物学

计算生物学

作者:互联网

计算生物学ex2

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

EXERCISE1

PART 1

在主目录下创建工作目录,用于存放此次作业所需要的文件

mkdir ex2

使用Xtpf7上传文档
上传之后展示如下
在这里插入图片描述

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文件,打开显示如图
在这里插入图片描述
表达基因可视化
包的安装同上
代码如下

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图在这里插入图片描述

标签:文件,计算,查看,ReadsPerGene,fastqc,tab,生物学,out
来源: https://blog.csdn.net/hy_628/article/details/111207292