author: "余顺太"
date: "2020-03-23"
knitr::opts_chunk$set(echo = TRUE)
主要参考的生信技能树文章：
原创10000+生信教程大神给你的RNA实战视频演练
上游流程一般不用学，公司测好序一般会直接给表达矩阵。
NCBI搜索文章 Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing
进入Data availability部分，拿到GSE的id号——点开超链接GSE111229进入GEO数据库，
dqvo7t.png
准备工作
- 安装conda
在中国大陆的小伙伴，需要更改镜像源配置
$ 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
- 安装软件
为了避免污染linux工作环境，推荐在conda中创建各个流程的安装环境，比如：
$ conda create -n rna python=2 #创建名为rna的软件安装环境
$ conda info --envs #查看当前conda环境
$ conda activate rna #激活conda的rna环境，曾健明使用的source activate rna，并不好用。
------------------以下为补充内容------------------
转录组分析需要用到的软件列表
质控：fastqc, multiqc, trimmomatic, cutadapt, trim-galore
比对：star, hisat2, bowtie2, tophat, bwa, subread
计数：htseq, bedtools, deeptools, salmon
------------------以上为补充内容------------------
conda install -y sra-tools multiqc trim-galore subread hisat2
$ conda search packagename #安装之前先检索是否存在该软件
$ 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
$ conda deactivate #注销当前的rna环境
转录组流程
第一步：sra数据下载和转fastq格式
下载SRA数据，去GEO里面找到SRA再找到文件。使用SRA Toolkit的prefetch进行下载，prefetch放在sratoolkit文件夹下的bin目录。
1.下载SRA Toolkit
SRA Toolkit的下载如下：
dqxqV1.png
找到对应的版本，右键复制链接地址，linux命令行使用wget下载
2.下载sra数据
prefetch下载数据，先将sratoolkit写入环境变量，然后直接下载
sra数据的下载浪费了老子一天的时间也没下载下来，去他马勒戈壁的，不下啦！(关掉电脑过了一夜，删掉了sratoolkit那个文件夹又TMD可以下了。)
$ conda activate rna
$ prefetch SRR6791442 #可以直接下载，下载的结果在~/ncbi/public/sra里面，这个文件夹会被自动创造出来。
也可以批量下载
$ cat srr.list
SRR6791441
SRR6791442
SRR6791443
SRR6791444
SRR6791445
$ cat srr.list |while read id; do (nohup prefetch $id &); done #曾健明说因为有管道命令，所以do后用()括起来。亲测有效，如果不写入bashrc则prefetch需要使用绝对变量。
$ ps -ef|grep prefetch #可以看到正在下载的东西
$ ps -l #上一步执行之后使用jobs -l没有反应，使用ps-l会有反应。 #事实上这两个都不太准确，应该用top命令。
$ less nohup.out|grep failed #虽然自动下载，但是有些是错的，需要重新下载。
我一个一个下的(这种方法下载的不是.sra结尾，但是文件是一样的)：
dqxl9O.png
每一个细胞都是一个单独的sra文件，是单细胞单样本；10X是只有一个样本，但是这个样本有4000到8000个细胞，但是这么多样本只有一个fastq文件，通过UMI和barcode去把细胞分开。
3.sra格式转fastq格式
单个sra格式文件转fastq格式
$ fastq-dump SRR6791442.sra #格式转还用到的软件是fastq-dump，该软件在sratoolkit的bin目录下。但是如果用conda装的软件，直接激活conda环境就可以直接用了。
sra格式批量转fastq格式
$ ls ./*sra |while read id; do (nohup fastq-dump --split-3 --gzip --skip-technical --clip ${id} &); done #曾健明说因为有管道命令，所以do后用()括起来。 #因为是单端测序所以只会转一个文件，若是双端测序会产生两个。
几个参数的意思(查一下~)：
--split-3：一个样本被拆分成两个fastq文件
--gzip：可以产生fastq的压缩文件，更省空间。
--skip-technical：
--clip：
-O #可以选择输出的位置，O是大写的。
曾健明说转了之后批量做比对，比对用hisat，所以要找到参考基因组的索引。如何下载参考基因组？？？
第二步：质量控制
1. fsatqc
对单个文件进行fastqc
$ fastqc SRR1039510_1.fastq.gz #对SRR1039510_1.fastq.gz进行fastqc。
$ fastqc -t 10 SRR1039510_1.fastq.gz #-t是线程，多线程更快一些。
$ fastqc -t 10 SRR1039510_1.fastq.gz -o /data1/lixianlong/YST/ #-o将结果输出到指定位置。
文件批量fastqc
$ ls *gz | xargs fastqc -t 2 #每个id_fastqc.html（比如SRR6791441_fastqc.html）都是一个质量报告
xargs命令可以将前面的文件作为参数传给后面的命令，也可以用while read id do，不可以用管道，因为管道传的是文件名不是文件。
2. multiqc
fastqc生成质控报告，multiqc将各个样本的质控报告整合为一个。
$ multiqc ./ #这一步生成的multiqc_report.html，是所有样本的整合报告，要确保已经安装了multiqc。
$ zless SRR1039508_1.fastq.gz #可以看到length是63，也就是算上adapter的长度是63。
在SRR1039508_1_fastqc.html文件里的Overrepresented sequences那一项里面可以找到adapter
dqzNM4.png
$ zcat SRR1039508_1.fastq.gz|grep GTCTTCTGCTTG
fastqc质控报告如何看自己去查！
第三步: 去除低质量reads和去除接头
运行如下代码，得到名为config的文件，包含两列数据(双端测序需要的)
$ 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
1.对单独一个文件进行trim-galore (如果trim-galore安装在了其他conda环境中，执行下列操作前需要先激活该环境。执行完之后再失活该环境）.
bin_trim_galore=trim_galore #这一步本是多此一举，但是很多时候并不想将软件写入环境变量，就可是在这里使用绝对路径，之后引用变量就可以了。
dir='/data1/lixianlong/YST/clean' #将结果导向的位置。
fq1='~/ncbi/public/sra/SRR1039508_1.fastq.gz'
fq2='~/ncbi/public/sra/SRR1039508_2.fastq.gz'
$bin_trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2
可以将得到的文件再进行fastqc和multiqc，和trim-galore之前的fastqc以及multiqc结果进行对比就可以看出差别，如果觉得不太满意可以再去调节参数。trim-galore之后得到的是干净的测序数据，就可以用来跑后面的流程了。
强烈警告：这个软件的名字叫做trim-galore，conda安装的使用也是使用 conda install trim-galore，但是安装成功之后可执行文件是trim_galore，并不是trim-galore，二者中间的横是不一样的，这点一定注意！否则无法执行。
--quality：设定Phred quality score阈值，默认为20。
--phred33：：选择-phred33或者-phred64，表示测序平台使用的Phred quality score。
--adapter：输入adapter序列。也可以不输入，Trim Galore! 会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
--stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
--length：设定输出reads长度阈值，小于设定值会被抛弃。因为前面发现算上adapter的长度是63，所以这里的length将阈值设为36，表示去掉两端adapter之后，长度小于36的会被丢弃。而如果算上adapter的长度是150，一般会将这里的length设为50，意思是将两端adapter去掉之后，长度小于50的会被丢弃。
-e：不是特别所谓
--stringency：adapter比对到多少算adapter，这里可以写成3，4，5都可以，后期可以慢慢调。
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
--output_dir：输入目录。需要提前建立目录，否则运行会报错。
-- trim-n : 移除read一端的reads
2.使用循环对文件进行批量的trim-galore操作
$ cat > qc.sh
bin_trim_galore=trim_galore #这一步本是多此一举，但是很多时候并不想将软件写入环境变量，就可是在这里使用绝对路径，之后引用变量就可以了。
dir='/data1/lixianlong/YST/clean' #将结果导向的位置。
cat $1 |while read id
do
arr=(${id}) # $id 表示$1的每一行，这里是config的每一行，以空格键分割。
fq1=${arr[0]} # 空格键分割的第一个元素
fq2=${arr[1]} # 空格键分割的第二个元素，曾健明说直接拿来用就好，为什么这样他也不知道。余顺太发现这个其实属于《LINUX命令行与shell脚本》6.7 数组变量的内容。
nohup $bin_trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
$ bash qc.sh config #config就是qc.sh脚本的$1
第四步：比对
为什么需要比对？
NGS测序下来的短序列（read）存储于FASTQ文件里面。虽然它们原本都来自于有序的基因组，但在经过DNA建库和测序之后，文件中不同read之间的前后顺序关系就已经全部丢失了。因此，FASTQ文件中紧挨着的两条read之间没有任何位置关系，它们都是随机来自于原本基因组中某个位置的短序列而已。因此，我们需要先把这一大堆的短序列捋顺，一个个去跟该物种的 参考基因组【注】比较，找到每一条read在参考基因组上的位置，然后按顺序排列好，这个过程就称为测序数据的比对。这也是核心流程真正意义上的第一步，只有完成了这个序列比对我们才有下一步的数据分析。
【注】参考基因组：指该物种的基因组序列，是已经组装成的完整基因组序列，常作为该物种的标准参照物，比如人类基因组参考序列，fasta格式。
序列比对本质上是一个寻找最大公共子字符串的过程。BLAST使用的是动态规划的算法来寻找这样的子串，但在面对巨量的短序列数据时，类似BLAST这样的软件实在太慢了！因此，需要更加有效的数据结构和相应的算法来完成这个搜索定位的任务，BWA就是其中最优秀的一个，它将BW(Burrows-Wheeler)压缩算法和后缀树相结合，能够让我们以较小的时间和空间代价，获得准确的序列比对结果。
为参考基因组构建索引——这其实是在为参考序列进行Burrows Wheeler变换（wiki: 块排序压缩），以便能够在序列比对的时候进行快速的搜索和定位。
$ bwa index human.fasta
以我们人类的参考基因组（3Gb长度）为例，这个构造过程需要消耗几个小时的时间（一般3个小时左右）。完成之后，你会看到类似如下几个以human.fasta为前缀的文件：
.
├── human.fasta.amb
├── human.fasta.ann
├── human.fasta.bwt
├── human.fasta.pac
└── human.fasta.sa
这些就是在比对时真正需要被用到的文件。这一步完成之后，我们就可以将read比对至参考基因组了：
转录组比对hisat2,subjunc,star这三个工具比较流行，基因组比对bwa,bowtie2比较流行；比对最重要的是用参考基因组构建索引，现有的参考基因组存储网站三个：ENSEMBL，UCSC，NCBI。每个软件构建索引的方式并不一样。现有比对工具在做mapping之前，都需要下载对应物种的参考基因组然后构建index。索引构建，自己去查一下各个软件构建索引的代码。索引构建非常慢，至少两个小时。
ENSEMBL的命名规则则是采用GRCh/m的方式，GRCh37对应hg19，hg38对应GRCh38。
1.参考基因组下载
1.打开USCS官网：http://genome.ucsc.edu/
Downloads —— Genome Data —— Human —— 找到Feb. 2009 (GRCh37/hg19)那里——点击Genome sequence files and select annotations (2bit, GTF, GC-content, etc) ——找到chromFa.tar.gz —— 右键复制链接地址wget下载
dqzyRO.png
下载hg19：（YST用的）
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz #可以新建一个文件夹叫 genome/hg19 ，然后下载在该文件夹中
tar -zxvf chromFa.tar.gz
cat *.fa > hg19.fa #对解压后的.fa文件进行合并, 这一步将所有的染色体信息整合在一起，重定向写入hg19.fa文件，得到参考基因组
rm chr*.fa #合并完成后删除合并前文件。
参考基因组注释下载(用于构建index)？？？
进入人和小鼠基因组注释信息官网GENCODE，选择data->human->GRCh37-mapped Releases，下载最新第26版本的hg19人类基因组注释信息。点击进入下载页面，将GTF和GFF3全部下载，解压。
其它的参考基因组下载方式
$ wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz #小鼠基因组，叫 genome/mm10，然后下载在该文件夹中.
$ wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz #可以新建一个文件夹叫 genome/hg38 ，然后下载在该文件夹中。
$ gunzip hg38.fa.gz # 使用gunzip进行解压。
2.建立索引
hisat索引文件下载（YST用的）
人类和小鼠的索引有现成的，HISAT2官网可以直接下载进行序列比对。HISAT2是TopHat2/Bowti2的继任者，使用改进的BWT算法，实现了更快的速度和更少的资源占用。
cd ~/reference
mkdir -p index/hisat && cd index/hisat
nohup wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz & # 我下的这个
nohup wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz &
nohup wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/grcm38.tar.gz &
nohup wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz & #小鼠的index
tar zxvf hg19.tar.gz
tar zxvf grcm38.tar.gz
tar zxvf hg38.tar.gz
有时候没有现成的index，我们就需要自己用HISAT2重新构建索引；包括外显子、剪切位点及SNP索引的建立。
hisat2构建索引方式自己去查
bowtie软件建立索引文件
cd ~/reference
mkdir -p index/bowtie && cd index/bowtie
nohup time ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ~/reference/genome/hg19/hg19.fa ~/reference/index/bowtie/hg19 1>hg19.bowtie_index.log 2>&1 &
nohup time ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ~/reference/genome/hg38/hg38.fa ~/reference/index/bowtie/hg38 1>hg38.bowtie_index.log 2>&1 &
nohup time ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ~/reference/genome/mm10/mm10.fa ~/reference/index/bowtie/mm10 1>mm10.bowtie_index.log 2>&1 &
bwa软件建立索引文件
cd ~/reference
mkdir -p index/bwa && cd index/bwa
nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index -a bwtsw -p ~/reference/index/bwa/hg19 ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1 &
nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index -a bwtsw -p ~/reference/index/bwa/hg38 ~/reference/genome/hg38/hg38.fa 1>hg38.bwa_index.log 2>&1 &
nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index -a bwtsw -p ~/reference/index/bwa/mm10 ~/reference/genome/mm10/mm10.fa 1>mm10.bwa_index.log 2>&1 &
3.mapping
3.1 单个进行的mapping
对每个文件的前1000行进行mapping(这里只是为了尝试一下，为了节省时间就只比对前1000行了。)
$ ls /data1/lixianlong/YST/clean/*.gz
$ ls /data1/lixianlong/YST/clean/*.gz|while read id; do (zcat $id |head -1000); done
$ ls /data1/lixianlong/YST/clean/*.gz|while read id; do (echo $(basename $id)); done|paste - - #与basename相对的是dirname，basename实现掐头功能，dirname实现去尾功能。
$ ls /data1/lixianlong/YST/clean/*.gz|while read id; do (zcat $id|head -1000 > $(basename $id)); done #将每个文件的前1000行留存下来，但是此时的名字虽然还是gz，其实并不是压缩文件了，可以直接用head和cat可以打开，这里只是没有改名字而已。文件名的后缀没有任何意义，真正起决定作用的是内容。可以新建一个文件夹运行这句代码，否则，源文件会被覆盖掉。
$ ls *|while read id; do (echo ${id%%.*});done #这样可以将.fq.gz都去掉，但是这里只是展示去掉之后的名字是什么，并没有真正去掉。
$ ls *|while read id; do (mv $id ${id%%.*});done #这里将所有文件都去掉名字后面的.fq.gz
$ ls *|while read id; do (mv $id $id.fq.gz);done #这里也可以将所有文件名字都加上.fq.gz
$ ls /data1/lixianlong/YST/clean/*.gz|while read id; do (zcat $id|head -1000 > $(basename $id ".gz")); done #也可以这样改掉，$id和".gz"一定要有空格，没有空格就是添加东西，有空格才是删掉。
对SRR1039508_1_val_1.fq.gz和SRR1039508_2_val_2.fq.gz的前1000行进行mapping
$ fq1=/data1/lixianlong/YST/clean/try/SRR1039508_1_val_1 #SRR1039508_1_val_1.fq.gz的前1000行保存在了文件SRR1039508_1_val_1里面
$ fq2=/data1/lixianlong/YST/clean/try/SRR1039508_2_val_2 #RR1039508_2_val_2.fq.gz的前1000行保存在了文件SRR1039508_2_val_2里面
$ hisat2_hg19_index='/data1/lixianlong/YST/reference/index/hisat/hg19/genome' #注意：hg19下面有很多文件，这里只写前缀genome！
$ hisat2 -p 10 -x $hisat2_hg19_index -1 $fq1 -2 $fq2 -S map.sam # 使用hisat2进行mapping，-S是将结果输出，mapping得到sam文件
$ hisat2 -p 10 -x $hisat2_hg19_index -1 $fq1 -2 $fq2 > map.sam # >和-S有同样功能。
使用hisat2比对得到的是97.40%的比对率，但是此处如果使用bwa比对率就是89.20%，比对率严重下降，因为hisat使用的是参考转录组，而bwa使用的是参考基因组，存在跨内含子的问题，导致比对率下降。事实上，对RNA-seq数据来说，不要使用bwa和bowtie这样的软件，它需要的是能进行跨越内含子比对的工具。
对SRR1039508_1_val_1.fq.gz和SRR1039508_2_val_2.fq.gz进行mapping
$ fq1=/data1/lixianlong/YST/clean/SRR1039508_1_val_1.fq.gz
$ fq2=/data1/lixianlong/YST/clean/SRR1039508_2_val_2.fq.gz
$ dir=/data1/lixianlong/YST/mapping_result
$ hisat2_hg19_index='/data1/lixianlong/YST/reference/index/hisat/hg19/genome' #注意：hg19下面有很多文件，这里只写前缀genome！
$ hisat2 -p 10 -x $hisat2_hg19_index -1 $fq1 -2 $fq2 -S $dir/SRR1039508.sam # 使用hisat2进行mapping，-S是将结果输出，使用>也可以，mapping得到sam文件
可以将以上代码包装进mapping.sh，然后qsub -cwd -l vf=8g,p=8 -V mapping.sh就可以了，qstat可以查询状态，qdel可以删除作业。
使用其他比对软件要用其它比对软件的index，index不是通用的。
$ ls *.sam|head #这个head的是管道前的output
$ ls *.sam|xargs head # 这个head的是前面输出结果的每一个变量所代表的文件，一定要掌握xargs的用法，有空搜索一下。
$ ls *.sam|xargs -i head {} # ls *.sam|xargs head 的完整写法
3.2 批量mapping
$ ls *gz|cut -d"_" -f1|sort -u|while read id; do
$ ls -lh ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz
$ hisat2 –p 10 –x /data1/lixianlong/YST/reference/index/hisat/hg19/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz –S ${id}_hisat.sam
$ done
4.sam格式转bam格式
批量转换
$ ls *.sam|while read id; do (samtools sort -O bam -@ 5 -o $(basename $id ".sam").bam $id);done #sam文件转成bam文件，文件变的更小，sort是很耗时间的。-@表示线程数量
$ top -u lixianlong # %CPU那一列就可以看出线程数，如果是五个线程，则状态最好的时候%CPU那一列值会接近500.
$ ls *bam|xargs -i samtools index {} #为bam文件建立索引。{}类似于$id，就是ls *bam的每一个结果会通过xargs传给后面的参数，这个参数就是{}。
#排好序并构建好索引，之后就可以使用IGV看了。
$ samtools view SRR1039508.bam |less -S # bam文件的查看使用samtools view，因为bam文件是按染色体sort过的，
$ samtools tview 08.bam #曾健明说这个用不着，平时都用IGV，因为IGV是可视化的，所以跳转很方便。
$ / #会出现goto，随便输入一个位置chr1:65891003看下，通过这种方式可快速到达指定位置。
$ samtools tview --reference /data1/lixianlong/YST/reference/genome/hg19/hg19.fa 08.bam #可以看到比对的情况，输入 chr1:220142425再试下。
$ samtools flagstat 08.bam # 也可以对bam文件进行统计，得到比对率之类的信息。samtools flagstat不能对sam进行统计，只能对bam统计。
sam格式怎样看自己去看微信公众号之类的里面的介绍。
转成bam之后将sam删除掉，因为sam太占内存，而且转成bam之后sam就没啥用了，bam本身就是sam的压缩。
曾健明说bam文件其实也可以进行qc
5.对bam文件再次fastqc
$ ls *bam|while read id; do (samtools flagstat -@ 10 $id > $(basename ${id} ".bam").flagstat ); done #实现批量操作，-@指定线程数量。
$ multiqc ./
第五步: counts
counts的原理就是比较基因上各个坐标对应的reads的个数多少。
# 如果一个个样本单独计数，输出多个文件使用代码是：
featureCounts -T 5 -p -t exon -g gene_id -a
# 如果一个个样本单独计数，输出多个文件使用代码是：
$ 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
$ done
# 如果是批量样本的bam进行计数，使用代码是：
$ 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 & # -T线程，-p双端测序，-t exon根据外显子来count，-g gene_id 表示count之后取基因名，得到的 all.id.txt文件就是表达矩阵，这个featureCounts有很多参数可调。
$ source deactivate
曾健明说对count也可以跑multiqc，就会生成一个报告，告诉你每个样本的比对率，比对上了多少兆reads。每一步都有qc，fastq文件有qc，bam文件也可以qc，count之后又有qc。表达矩阵拿到R里面去做差异分析也有qc，每一步都要搞清楚参数有没有用对，qc无处不在！
有了这个表达矩阵之后所有的下游分析就可以开始了。
第六步：DEG
1. 对得到的表达矩阵进行初步探索
rm(list = ls()) #清空环境变量，读表达矩阵要做的第一步。
options(stringsAsFactors = F) #读表达矩阵要做的第二步。
a=read.table('all.id.txt',header = T) #如果不加header = T，则表头位置是V1，V2，V3...，为了把列名放在表头的位置要加header = T。
tmp=a[1:14,1:7]
meta=a[,1:6] #左边的6列是meta信息，meta信息就存着每一个基因的名字，坐标和长度。
exprSet=a[,7:ncol(a)] #7列之后的是表达矩阵信息，
colnames(exprSet)
a2=exprSet[,'SRR1039516.hisat.bam']
#因为曾健明使用bowtie，bwa，hisat，subjunc这几种比对软件，他做了一下相关性分析，看一下几种比对软件结果的相关性。
library(corrplot) #看相关性的第一种方法。
png('cor.png') #打开一个画板
corrplot(cor(log2(exprSet+1))) #log多少其实没有太所谓，+1的目的是为了防止0值。 #作图
dev.off() #关闭画板。
library(pheatmap) #看相关性的第二种方法。画heatmap的包有十多个，这里只用了这一个。
png('heatmap.png')
m=cor(log2(exprSet+1))
pheatmap(scale(cor(log2(exprSet+1)))) #进行一下归一化，这样热图差异更明显。曾健明的结果发现同一个软件的比对结果聚集在了一起，说明这些软件的差异性要比样本之间的差异性要大。
dev.off()
曾健明说：转录组上游基本上就是在讲linux命令，转录组下游就是在讲R语言。转录组和外显子组测序分析本身就不存在，只不过是利用linux和R的命令来分析数据。
2. 使用airway数据包做示范
library(airway) #airway是个数据包，从Bioconductor上面下载，数据包就是整个包就是为了提供一个数据。
data(airway)
airway # 我们想取airway的表达矩阵，可以先看一下airway，因为airway是一个表达矩阵。
?RangedSummarizedExperiment
exprSet=assay(airway) #发现使用assay可以取到表达矩阵。
colnames(exprSet)
# 首先对这个表达矩阵看一下相关性
library(corrplot) #看相关性的第一种方法。
png('cor.png') #打开一个画板
corrplot(cor(log2(exprSet+1))) #log多少其实没有太所谓，+1的目的是为了防止0值。 #作图
dev.off() #关闭画板。从右下角Files里面可以看到图,发现相关性都接近于1
library(pheatmap) #看相关性的第二种方法。画heatmap的包有十多个，这里只用了这一个。
png('heatmap.png')
m=cor(log2(exprSet+1))
pheatmap(scale(cor(log2(exprSet+1)))) #进行一下归一化，这样热图差异更明显。曾健明的结果发现同一个软件的比对结果聚集在了一起，说明这些软件的差异性要比样本之间的差异性要大。
dev.off()
#是否符合样本的实际特征呢？
tmp=colData(airway)
tmp=as.data.frame(tmp) #dex那一列点一下可以排序，发现SRR1039509，SRR1039513，SRR1039517，SRR1039521这四个是处理，另外四个是不处理的，从上面pheatmap的结果可以看到SRR1039509，SRR1039513,SRR1039521,但是SRR1039517却没有聚在一起。
# 也可以使用hclust去查看是否符合样本的实际特征
## hclust
group_list=colData(airway)[,3]
colnames(exprSet)=paste(group_list,1:ncol(exprSet),sep='_')
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), cex = 0.7, col = "blue")
hc=hclust(dist(t(log2(exprSet+1))))
par(mar=c(5,5,5,10))
png('hclust.png',res=120)
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
dev.off() #从图上的结果可以看到trt_6没有聚到一起，trt_6就是SRR1039517.如何确定trt_6就是SRR1039517? 通过view(tmp)就可以看到原始的排序，trt6对应的就是SRR1039517.从avgLength那一列也可以看出来，SRR1039517是87，而另外三个trt都是126
exprSet=assay(airway) #使用assay可以取到表达矩阵。
colnames(exprSet)
a1=exprSet[,'SRR1039516']
a1 # a1就是每个基因的counts数
group_list=colData(airway)[,3]
3. 使用airway数据包产生的a1和all.id.txt(我们自己做的矩阵)产生的a2进行比较。
a2=data.frame(id=meta[,1],a2=a2) # a2是自己做出来的数据
a1=data.frame(id=names(a1),a1=as.numeric(a1)) # a1是airway数据包产生的，a1=as.numeric(a1)表示只要它的数字内容。
library(stringr) #使用stringr包的str_split可以实现分割。
a2$id <- str_split(a2$id,'\\.',simplify = T)[,1] # 因为a2的基因名字是有问题的，比如ENSG00000237613.2,所以要将.2去除，只留前面的名字。
# 对a2的id列进行分割，\\.表示以点为分割位置，[,1]表示分割好之后取分割好的第一列。
tmp=merge(a1,a2,by='id') #可以根据id将a1和a2合并
png('tmp.png')
plot(tmp[,2:3])
dev.off() #从得到的图上可以看到有一些离群点的存在。
一定要看的文章(晚上服务器比较好用，先不看了，明天看。2020年3月27日)
从零开始完整学习全基因组测序（WGS）数据分析：第1节 DNA测序技术
从零开始完整学习全基因组测序（WGS）数据分析：第2节 FASTA和FASTQ
从零开始完整学习全基因组测序（WGS）数据分析：第3节 数据质控
从零开始完整学习全基因组测序（WGS）数据分析：第4节 构建WGS主流程
网友评论