control+左右键：光标在单词之间跳转
control+a：光标跳到首
control+e：光标跳到尾
more 查看文件内容，且显示前几行。q 退出，文本可见
less q退出，文本不可见
cat 打开所有内容，文件太大了不可用。
du -h NAME 查看文件大小（人类看得懂的）
ls ..
cd ..：..表示返回上一层目录。一个.表示当前目录
cd -：返回刚才的目录
gzip：压缩文件
ls *bam:查看bam文件
’>>‘ 不会覆盖原来的文件，在原来的文件后面添加东西
‘>’ 覆盖原来的文件
sra——fastq——fastqc
下载玉米参考基因组文件Fastq
Ensemble网站上Zea Mays的参考基因组文件下载，注意要选择toplevel，染色体数目全的，而不是单条染色体的文进组文件。
下载玉米基因注释文件GFF
同上
MAC系统：ensemble网站下载参考基因组，显示在finda里，找到目标文件可以先解压，通过FileZilla上传到服务器
FileZilla用法：
连接服务器，输入ip地址，输入qwang，输入密码，端口21.快速连接——找到MAC里的文件——选择linux里的要保存的文件夹——双击MAC里的文件自动上传
亚硫酸盐测序=甲基化测序
数据哪里找
文章的DATA部分，数据存放在NCBI的GEO数据库或SR数据库
SRP号包含所有数据
打开NCBI——找到GEO数据库——搜索GEO号——找到RNAseq的数据点击GSE号——samples里有5个，3个生物学重复，两个对照——点击其中一个重复——找到SRX号——找到SRP号——点击RUN：5——全选5个SRR号文件数据——Download选择Metadata——下载一个csv文件——这个table文件涵盖
已知SRR号——百度搜索如何在NCBI网站上下载数据——找到FTP下载的路径——修改末尾SRR号——输入wget - r ftp连接(这个方法我没成功)
FTP
一个数据库，MAC OS系统下做成OS版的可视化，显示为finda。FTP里有各种文件夹，SRR/SRR191/SRR1916145/意思是SRR文件夹里SRR191文件夹里的SRR1916145这个文件夹，里面想要的文件是SRR1916145.sra
sra toolkit的下载
https://github.com/ncbi/sra-tools/wiki/HowTo:-Binary-Installation
找到Ubuntu，右键复制链接，wget+链接下载
bin（binary的简写）
下载一个.tar.gz的压缩文件，用命令tar zxvf +sratoolkit.tar.gz解压文件，解压之后的文件里有bin文件，打开bin文件发现有很多命令，我们用fastq-dump命令
在Mac里下载sratoolkit同理。（Mac里的sratoolkit我还不会用）
（NCBI下载SRA数据的三种方法：FTP，迅雷，sratoolkit）
(1)用迅雷
屏幕快照 2020-05-06 下午8.36.32.png
屏幕快照 2020-05-06 下午8.41.17.png
复制链接，复制到迅雷用迅雷下载，或者wget下载。但是考虑到网速，不知道有没有更快的方法。
下载好的文件通过FileZilla传输到服务器上
(2)用sratoolkit
将sratoolkit安装到服务器，用prefetch的命令
prefetch SRR19116154
对文件进行格式转化
~/biosoft/sratoolkit.2.10.5-ubuntu64/bin/fastq-dump SRR1916154.1
将sra格式转化成fastq格式
屏幕快照 2020-05-06 下午10.13.48.png
less 查看文件，q退出
可以将prefetch 和fastq-dump命令写入环境变量
修改环境变量
vim ~/.bashrc
下拉到最后，在文件最后一行写入：
export PATH=命令的绝对路径/bin:$PATH
esc退出，:wq!保存
最后激活一下修改后的bashrc
source ~/.bashrc
echo $PATH查看所有环境变量
将一个文件夹从一个服务器传输到另一个服务器
scp -r #文件夹路径# 目标用户名@ip地址：home
注意：tab补齐路径，在正确的路径下才会有效，如果tab失灵说明路径不对。同理正确的命令tab才会补齐，如果失败那么说明命令不对。
质控
将fastqc添加到环境变量,直接进输入命令：
fastqc SRR1916154.fastq
得到一个SRR1916154.fastqc.zip
解压zip文件的命令：
unzip SRR1916154.fastqc.zip
得到一个HTML格式的质检报告（QC report）用网页打开看
因为酵母这个数据质检结果质量还不错，直接用来比对
建立index，以便序列比对
1）下载一个hisat2
bing搜索hisat2 download
找到linux版本，右键复制连接，wget下载
或者复制连接，迅雷下载到本地，本地上传到服务器，unzip解压。
将hisat2安装在biosoft文件夹里，ls --color查看hisat2是蓝紫色的，可以应用。
2）将hisat2写进环境变量。
3）建立index
用hisat2-build命令
hisat2-build Saccharomyces_cerevisiae.R64-1-1.dna_sm.toplevel.fa#酵母的参考基因组# yeast_ref#输出的文件的名字#
hisat2 --help命令查看命令的用法，<>里是必填项
因为酵母基因组很小，hisat2-build建立ind很快，如图：
屏幕快照 2020-05-11 下午4.03.15.png
输出名为yeast_ref的文件8个
4）建立一个文档，记录建立索引的这个命令
vi cmd1_build_index.sh
点击i，此时是输入文本的状态，输入hisat2-build Saccharomyces_cerevisiae.R64-1-1.dna_sm.toplevel.fa yeast_ref记录本次命令，下次忘记或者怎样可以用来查看，查看的命令是more
5）将实验数据同index进行比对
（1）mkdir analysis/read_aln新建一个储存分析结果的文件夹，在此文件夹下hisat2
（2）文件SRR1916154.1.fastq所在路径：~/data/yeastproject/data/RNAseqDATA
（3）index所在路径：~/data/yeastsqence
准备好之后，执行命令：
hisat2 -x ~/data/yeastsqence/yeast_ref#index的路径# -U ~/data/yeastproject/data/RNAseqDATA/SRR1916154.1.fastq#质控后的实验数据# >
上述命令是标准结构的命令，但将结果直接打到屏幕上，并不是我们想要的。
修改后为以下命令：
nohup hisat2 -x ~/data/yeastsqence/yeast_ref -U ~/data/yeastproject/data/RNAseqDATA/SRR1916154.1.fastq > 3b_strain_2.sam &
nohup#挂在后台运行# time#显示运行时间# hisat2 -x ~/data/yeastsqence/yeast_ref -U ~/data/yeastproject/data/RNAseqDATA/SRR1916154.1.fastq > 3b_strain_2.sam#输出文件命名为这个# &#与nohup构成完整命令#
注意：time在基因组庞大的数据下不能用，不然后期改很麻烦。因为这个软件输出的分析后的文件格式是sam，直接命名.sam
（4）在analysis/read_aln/文件夹下输入jobs查看正在运行的文件
输入top -c显示运行的状态
输入tail 3b_strain_2.sam查看写到哪里了，tail -1 3b_strain_2.sam查看书写的最后一行，根据最后一行不断更新变化判断是否正在进行比对。jobs命令显示DONE，表示完成：
屏幕快照 2020-05-11 下午5.54.58.png
注意：有一个问题，当我们用time命令的时候，会把时间写进去，（写在了最后一行，tail可看）导致文件格式出现错误。所以，不要time这个命令重新运行一遍（因为视频教程是如此操作，现在看来这步浪费时间，运行基因组庞大的可不能使用time命令），重新运行后覆盖原文件。
（5）完成后，查看3b_strain_2.sam
less -s 3b_strain_2.sam
sam文件比较占地方，用samtools软件处理一下。
samtools软件的安装：conda install samtools
用到samtools view 命令
samtools view -bS 3b_strain_2.sam -o 3b_strain_2.bam
将sam文件转化成bam文件。
或者在命令首位加上nohup &，在后台运行。
生成bam文件，按照染色体号来给文件内容排序，用samtools sort命令
samtools sort 3b_strain_2.bam -o 3b_strain_2.srt.bam
计数
（1）安装htseq软件
conda install htseq
（2）htseq-count ，下载酵母基因注释文件gff格式,或者gtf格式，只要把FTP地址中改成gff，gtf即可。
路径：~/data/yeastproject/ref/Saccharomyces_cerevisiae.R64-1-1.47.gff3
（3）新建一个文件夹read_count，在analysis目录下
(4)htseq-count
路径：~/analysis/read_aln/3b_strain_2.srt.bam
htseq-count -f bam -r pos ~/analysis/read_aln/3b_strain_2.srt.bam ~/data/yeastproject/ref/Saccharomyces_cerevisiae.R64-1-1.47.gtf > 3b_strain_2.count.tab
输出一个3b_strain_2.count.tab文件
将两个空白对照和三个样本数据都进行一个计数处理
（5）将5个计数后的count.tab文件进行合并：
用python写脚本的方式
在read_counts目录下新建一个scripts文件夹，在此文件夹里写入一个脚本文件
vi merge_read_count.py
建好这个文本文档，准备在里面输入编程程序
（6）read_count下，在linux命令行里输入
python scripts/merge_read_count.py EV_strain_3.count.tab,EV_strain_4.count.tab,3b_strain_2.count.tab,3b_strain_3.count.tab,3b_strain_4.count.tab EV_strain_3,EV_strain_4,3b_strain_2,3b_strain_3,3b_strain_4
意思是：python这个脚本文件（.py），把五个.tab文件依次对应输出名字为这几个名字。一一对应不要搞错。
每次修改.py里的内容，python出来的结果都会不一样。
（7）写程序
vi xxxxxxxxxx.py
import sys
sample_counts = sys.argv[1]
sample_names = sys.argv[2]
count_files = sample_counts.split(",")
sample_ids = sample_names.split(",")
#print(count_files)
count_dict = {}
for count_file in count_files:
with open(count_file) as count:
for line in count:
if line.startwith("__"):
continue
line = line.rstrip("\n")
ele = line.split("\t")
#print(ele)
if ele[0] in count_dict:
count_dict[ele[0]].append(ele[1])
else:
count_dict[ele[0]] = [ele[1]]
print("gene_id\t" + "\t".join(sample_ids))
for gene_id in count_dict:
#print(count_dict[gene_id])
print(gene_id + "\t" + "\t".join(count_dict[gene_id]))
try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
熟练掌握可用管道命令一步到位
hisat2 -x (yeast_ref,indexd的路径) -U （fastq原始文件路径）|samtools view -bS - |samtools sort - -o 3b_strain_2.srt.bam
python编程
下载数据数量庞大的时候，比如要下载一百个数据，需要编写脚本scripts
编程思路：教我老弟学生信第5天
RNAseq数据的IGV展示
（1）下载igv：bing搜索igv，第一个条点击download，找到linux版本，右键复制链接，wget下载，unzip压缩文件，cd 打开解压文件，ls查看里面的项目，sh igv.sh命令打开igv，弹出窗口。
（2）针对Mac，下载igv for Mac ，因为云服务器太小了。linux打不开
（3）在igv界面里，选择参考基因组，如果没有，通过上传的方式，上传toplevel文件。点击Genomes——from file选择文件。（在哪里安装的就从哪里上传）
（4）igv界面，点击file——load from file选择gtf文件。
（5）igv界面，点击file——loadfromfile选择bam文件。
（6）这样就可以看到可视化的比对结果，显示比对上的每一个reads。
（7）保存：file——save session
R处理
setwd("D:\\yeast")
count_tab <- read.table("yeast_EV_DNMT3B_count.tab", header = T)
rownames(count_tab) = count_tab$gene_id
count_tab = count_tab[, -c(1)]
colData <- read.table("sample_info4DESeq2.txt", header = T)
colData$condition = factor(colData$condition, c("EV", "DNMT3B"))
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_tab,
colData = colData,
design= ~ condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_DNMT3B_vs_EV")
res <- res[order(res$padj),]
resDF = as.data.frame(res)
write.csv(resDF, file = "yeast_DESeq2_DEG.csv")#把文件写成csv文件保存
（注释# or to shrink log fold changes association with condition:
res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")）
PCA图
library(ggplot2)
## PCA
vsd <- vst(dds, blind=FALSE)
pcaData <- plotPCA(vsd, intgroup=c("condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition)) +
geom_point(size=3) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed()
Rplot.png
两个颜色离得很远，分得很开，说明数据好，有差异
MA-plot
plotMA(res, ylim=c(-2,2))
MA-plot.png
远离红线的是差异表达基因，离得越远，差异越大
网友评论