#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文件 复制到电脑桌面看
