生信学习笔记
linux部分功能
查看文件夹
工具 选项 可以设置鼠标功能 可以设置右键粘贴
双击这个窗口可以再打开一个窗口 互不干扰
Rnaseq数据挖掘
观看陈魏学视频 了解RNAseq
测序技术介绍及基本流程
我们能做什么：火山图 热图差异基因分析
蛋白互作图
kegg
GO 生物学功能和表达位置
泡泡图 分析结构差异
这几个图回头还得再好好研究研究
学习材料：简书原创10000+生信教程大神给你的RNA实战视频演练
实例分析
安装conda
https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/ 下载需要的版本（conda2 linux版本）
更改镜像 翻墙也行应该
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
安装软件
创建一个软件安装环境 放置污染大环境
conda create -n rna python=2 # 环境名字是rna 可以自己随便改别的
查看conda环境 conda info --envs
工作时候载入环境source activate rna
安装软件
质控
fastqc 查看fastq文件质量
multiqc 合成不同样本的质控报告 一起查看
trim-galore 过滤低质量数据
trimmomatic, cutadapt 好像没用到
比对
star, hisat2（这个最好）, bowtie2, tophat, bwa, subread
计数
htseq, bedtools, deeptools, salmon
https://bioconda.github.io/recipes.html 这个地址可以查看conda是否有你要安装的软件
conda安装成功后 以下是安装RNAseq用到的软件的代码
conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
转录组流程
1. 下载sra数据
NCBI可以下载SRR list 里面有SRR号码 一个一排 也可以新建一个SRR_Acc_List.txt文档 也是一个号码一行。这个list是让代码读取SRR号码以在网上下载用的
下载SRA数据代码
wkd=/home/zhaowei/project/airway/ #设置工作目录
source activate rna 载入工作环境
cat SRR_Acc_List.txt | while read id; do (prefetch ${id} &);done
ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done #没太看懂这个 prefetch是下载SRA用的
2.将SRA格式转成fastq格式
用fastq-dump软件
单个文件转换
nohup fastq-dump --split-3 --skip-technical --clip --gzip $i & #把$i换成文件名
批量转换（循环）
for i in $wkd/*sra
do
echo $i
nohup fastq-dump --split-3 --skip-technical --clip --gzip $i &
done # 循环代码没看懂 nohup可以让进程在后台运行 不影响前台操作
因为现在大部分原始数据为双端测序 即从3和5两端测序 所以一个SRA文件可以转换成两个fastq文件即＿１和＿２，--split-3 这个参数就是转换双端测序文件用的 如果是单端测序文件 就不加这参数
生成的文件格式例子：SRR1039522_1.fastq.gz SRR1039522_2.fastq.gz
3.数据的质量检查
ls *gz | xargs fastqc -t 10 # 用fastqc进行检查 并且生成一个html的报告 用打开文件夹的那个按钮打开查看
multiqc ./ #整合生成的质控报告 html结果查看方法同上
4.过滤低质量数据
mkdir $wkd/clean
cd $wkd/clean
ls /home/jmzeng/project/airway/raw/*_1.fastq.gz >1
ls /home/jmzeng/project/airway/raw/*_2.fastq.gz >2
paste 1 2 > config
#把两个fastq.gz文件赋值 粘贴到config中并形成2列
用trim_galore过滤
source activate rna
bin_trim_galore=trim_galore #把bin_trim_galore赋值 下面有代码
dir='/home/zhaowei/project/airway/clean' #输出文件夹 自己设定
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]} # 如果fq1在第二列 把0换成1
fq2=${arr[1]} #道理同上
nohup $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
source deactivate #退出环境 也可以不退
#以上代码可以保存成 .sh结尾的文件 下次用bash一下就行 不用总写了 vim进去改改就行
5.比对
实际用时候感觉用一个软件就行 hisat2
1.运行单个样本 测试
mkdir $wkd/test
cd $wkd/test
source activate rna
ls $wkd/clean/*gz |while read id;do (zcat ${id}|head -1000> $(basename ${id} ".gz"));done
id=SRR1039508 # id写成需要的id
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}_1_val_1.fq -2 ${id}_2_val_2.fq -S ${id}.hisat.sam #注意看软件说明书 输出的是sam还是bam 这里写的只是我们自己起的名.sam 线程分配的数字可以改（-p 10那地方） 如果只用一个软件比对下面就不用了
subjunc -T 5 -i /public/reference/index/subread/hg38 -r ${id}_1_val_1.fq -R ${id}_2_val_2.fq -o ${id}.subjunc.sam
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38 -1 ${id}_1_val_1.fq -2 ${id}_2_val_2.fq -S ${id}.bowtie.sam
bwa mem -t 5 -M /public/reference/index/bwa/hg38 ${id}_1_val_1.fq ${id}_2_val_2.fq > ${id}.bwa.sam
2.批量代码
cd $wkd/clean
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat.sam
subjunc -T 5 -i /public/reference/index/subread/hg38 -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fq.gz -o ${id}.subjunc.sam
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38 -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.bowtie.sam
bwa mem -t 5 -M /public/reference/index/bwa/hg38 ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz > ${id}.bwa.sam
done
6.sam文件转换成bam
ls *.sam|while read id ;do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam ${id});done
rm *.sam
#报错时候看看sam文件到底转完没有 我这次就是比对到一半就运行了这个命令结果失败
7.建立bam文件索引
ls *.bam |xargs -i samtools index {} #每个文件都有索引哟
8.reads比对情况统计
ls *.bam |xargs -i samtools flagstat -@ 10 {} >
ls *.bam |while read id ;do ( nohup samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat & );done
source deactivate
最终结果示例：
9.计数
mkdir $wkd/align
cd $wkd/align
source activate rna #建文件夹 激活环境 （好像是可选）
一个一个样本计数 输出多个文件
for fn in {508..523}
do
featureCounts -T 5 -p -t exon -g gene_id -a /public/reference/gtf/gencode/gencode.v25.annotation.gtf.gz -o $fn.counts.txt SRR1039$fn.bam #SRR这地方似乎可以根据实际情况改一下
done
一起计数 输出一个文件
mkdir $wkd/align
cd $wkd/align
source activate rna
gtf="/public/reference/gtf/gencode/gencode.v25.annotation.gtf.gz"
featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o all.id.txt *.bam 1>counts.id.log 2>&1 &
# 代码看不懂 输出的all.id.txt文件就是表达矩阵 可以用来进行下一步差异分析画图等等 用R
count后文件示例
10.数据检查
rm(list = ls())options(stringsAsFactors =F)a=read.table('all.id.txt',header =T)tmp=a[1:14,1:7]meta=a[,1:6]exprSet=a[,7:ncol(a)]colnames(exprSet)a2=exprSet[,'SRR1039516.hisat.bam']library(airway)data(airway)exprSet=assay(airway)colnames(exprSet)a1=exprSet[,'SRR1039516']group_list=colData(airway)[,3]a2=data.frame(id=meta[,1],a2=a2)a1=data.frame(id=names(a1),a1=as.numeric(a1))library(stringr)a2$id <- str_split(a2$id,'\\.',simplify =T)[,1]tmp=merge(a1,a2,by='id')png('tmp.png')plot(tmp[,2:3])dev.off()library(corrplot)png('cor.png')corrplot(cor(log2(exprSet+1)))dev.off()library(pheatmap)png('heatmap.png')m=cor(log2(exprSet+1))pheatmap(scale(cor(log2(exprSet+1))))dev.off() #一点也看不懂 SRR应该是能改的
小TIPS
软链接 ln -s /路径/文件 ./ 这个./别忘了打了 表示链接到当前文件夹
vim后点i 可以编辑 退出时候 用shift：qw 退出并保存
每次用软件命令之前 看看自己文件是不是那么回事 less -SN 整齐的看 别用cat 太大了
设置好变量用echo $看一看 对不对
vim用 ：set paste命令可以用粘贴模式 格式无损
