#RNA-seq实战代码
#正常情况下公司给我们的数据就直接是fastq格式文件
#使用fastqc查看测序数据质量,就会得到*.html文件
fastqc *.fastq.gz
#在html文件所在目录下将*.html文件拷入一个单独的文件夹:cp *.html 新建的单独的文件夹所在的目录
#将linux中的文件复制到桌面
cp -r 要复制的文件所在目录 /mnt/c/Users/happy/Desktop/
#将所有的质量报告合并:对html文件所在目录跑一下multiqc:multiqc 新建的单独的文件夹所在的目录
multiqc 新建的单独的文件夹所在的目录 (要在rnaseq3.7环境中跑这个代码)
#使用trim_galore 过滤fastq文件中质量不好的reads和remove接头序列
ls /home/happy/project/airway/raw/*_1.fastq.gz >1
ls /home/happy/project/airway/raw/*_2.fastq.gz >2
#生成只包含*—1.fastq.gz和*—1.fastq.gz文件的目录
paste 1 2 >config
#过滤文件打算存储的位置
dir='/home/happy/project/airway/clean/'
#循环如下:
cat config |while read id;
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
nohup trim_galore -q 25 --phred33 --stringency 3 --length 36 -e 0.1 --paired $fq1 $fq2 --gzip -0 $dir &
done
#3、再次使用fastqc查看过滤后的文件
#比对,转录组比对主要是使用hisat2、subjunc、star;基因组主要是使用bwa和bowtie2
#比对首先要下载参考基因组构建索引index,也可以直接在官网下载index,人类和小鼠的index有现成的
#方法一:直接在hisat2官网下载index,选择合适的物种
#hisat2官网:http://daehwankimlab.github.io/hisat2/download/
#构建存储索引的文件夹
mkdir index
cd index
wget 复制链接
#解压index文件
tar -zxvf *.tar.gz
#删除压缩包
rm -rf *.tar.gz
#方法二:自己下载参考基因组和注释文件构建index
#下载基因组文件
wget 参考基因组链接
tar -zxvf *.tar.gz
#下载注释文件
wget 注释文件链接
gunzip *.gtf.gz
extract_exons.py hg19.ncbiRefSeq.gtf > genome.exon 其中hg19.ncbiRefSeq.gtf是解压后的注释文件
extract_splice_sites.py hg19.ncbiRefSeq.gtf > genome.ss
#建立index
hisat2-build hg19.fa --ss genome.ss --exon genome.exon genome_tran
#4.比对
hisat2 --dta -t -x /data/RNAseq/mm10/genome -1 /data/RNAseq/cleandata/trim_galoredata/CK-4_1_val_1.fq.gz -2 /data/RNAseq/cleandata/trim_galoredata/CK-4_2_val_2.fq.gz -S cleandata/hisat2_mm10data/CK4.sam
#/data/RNAseq/mm10/ 是参考基因组文件所在目录,genome这个不是文件夹,是index文件的前缀,mm10文件下并没有这个文件,如果不加genome就会发生报错
#-1和-2分别表示双端测序的1个文件,后面跟的是文件路径
#-S 后面跟的是输出数据所在文件夹
#-x 后面接的是索引文件
#命令没有问题的话,出现下面提示表示程序正在执行,就去干其他的,等执行结束再往下:
Time loading forward index: 00:02:51
Time loading reference: 00:00:16
#写循环批量处理(for in 循环和while循环各写一个)
#for in 循环,创建脚本文件
vim hisat2_mm10_batch.sh
#输入下面内容:
for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome
-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz
-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz
-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam;done
#保存脚本,再运行脚本
bash hisat2_mm10_batch.sh
#while 循环
#先创建存储了SRR号的文件。例如sra.txt
less sra.txt|while read id ;
do
hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome
-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz
-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz
-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam;done
#subjunc比对会直接输出bam文件,就不用后续的sam转bam文件
subjunc -T 5 -i /home/happy/project/airway/reference/mm10/genome -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fa.gz -o ${id}.subjunc.bam
#5.samtools将sam文件转bam文件,samtools可以进行格式转换、排序、索引
#第一步:转换
#SRR.txt为存储了SRR号的文件
单样本处理:samtools view -@8 -b SRR1039509.sam > SRR1039509.bam
多样本处理:less /home/happy/project/airway/hisat2/SRR.txt|while read id ;do samtools view -@8 -b {$id}.sam > {$id}.bam done
#第二步:排序 samtools sort
#其中 -l是(-L的小写)
单样本处理:samtool sort -@ 8 -l -o SRR1039508.bam.sort SRR1039508.bam
多样本批量处理:ls /home/happy/project/airway/SRR.txt |while read id ;do samtools sort -@ 8 -l {$id}.bam -o {$id}.bam.sort;done
#第三步:建立索引(index命令)在SRR1039508.bam.sort所在文件夹下操作
单样本处理:samtools index SRR1039508.bam.sort SRR1039508.bam.index
查看sam文件直接用less,查看bam文件用samtools view SRR1039508.bam |less -S(分页查看) 可以看出来是按pos排序还是按name排序,这里是按pos(染色体位置)
#第四步:查看reads比对情况(flagstat命令),输出.falgstat文件
samtools flagstat -@ 8 SRR1039508.bam.sort > SRR1039508.bam.flagstat
#第五步:igv可视化比对结果
先将bam文件转为bw文件,要先构建索引 samtools index SRR1039508.bam.sort
再使用bamCoverage转换格式:bamCoverage -b SRR1039508.bam.sort -o SRR1039508.sort.bw
再将bw文件复制到桌面:cp *.bw /mnt/c/Users/happy/Desktop/ (bw文件要指定文件所在目录)
igv里面选择file 载入bw文件
#6.计数(htseq、featureCounts )通常用featureCounts计数
htseq计数代码(很慢)
htseq-count -f bam -r pos /home/happy/project/airway/hisat2data/SRR1039508.bam.sort /home/happy/project/airway/reference/gencode.v10.annotation.gtf.gz
#-f 指定bam文件格式
#-r 指定排序方式,这里是按位置position
#/home/happy/project/airway/hisat2data/SRR1039508.bam.sort 进行计数的bam文件
#/home/happy/project/airway/reference/gencode.v10.annotation.gtf.gz 注释文件所在目录
featureCounts计数(很快):参考基因组和GTF/GFF/SAF注释文件来自同一个网站,同一个版本
featureCounts 需要两个输入文件:比对产生的BAM/ SAM文件(一般会用bam文件,因为所占空间小);区间注释文件(GTF格式, SAF格式)
代码1:featureCounts -T 10 -a ../reference/gencode.v10.annotation.gtf.gz -o count.txt -p -B -C -f -t exon -g gene_id /home/happy/project/airway/hisat2data/*sort
代码2:featureCounts -T 10 -a ../reference/gencode.v10.annotation.gtf.gz -o count1.txt -p -B -C -t exon -g gene_name /home/happy/project/airway/hisat2data/*sort
-g 后面也可以接gene_name
#如果代码1比对率不高,说明不是全外显子测序,则用代码2
-T 指定线程数 -f 如果-f被设置,那将会统计feature层面的数据,如exon-level,否则会统计meta-feature层面的数据,如gene-levels(feature数目指的是exon数目 meta-feature 指的是基因数目)
-a 参考GTF文件名 -F 参考文件的格式,一般为GTF/SAF ,默认为GTF -p 双端测序 -B 在-p选择的条件下,只有两端reads都被比对上的fragment才会被统计
-t 设置feature-type,-t指定的必须是gtf中有的feature,同时read只有落到这些feature上才会被统计到,默认是“exon”
-g 默认为gene_id -o 输出文件的名字,可输出被统计数据的txt文本及summary文本
/home/happy/project/airway/hisat2data/*sort 被统计的文件所在目录及名称
#featureCounts查看及结果解析:less -SN count.txt (-SN可以使查看是排版更好)
Geneid:基因的ensemble基因号; Chr:多个feature所在的染色体编号; Start:多个feature起始位点,与前面一一对应; End:多个feature终止位点,与前面一一对应
Strand:正负链 Length:基因长度 sampleID:一列代表一个样本,数值表示比对到该基因上的read数目
feature数目指的是exon数目 meta-feature 指的是基因数目
#7.表达矩阵探索
multiqc 要在3.7版本的python才不会报错,所以要创建pyhton3.7环境下下载multiqc :conda create --name rnaseq3.7 python=3.7
conda activate rnaseq3.7
conda install -c bioconda -c conda-forge multiqc
multiqc count1.txt.summary 就会得到一个html文件 复制到电脑桌面看