软件安装安装bioconda:
https:///miniconda.html 
# 安装
bash Miniconda3-latest-MacOSX-x86_64.sh
之后会要求看一个license，和安装环境路径，输入yes 就好。 source .bash_profile
# 查看是否安装成功
conda
出现上图代表安装成功 conda config --add channels conda-forge
conda config --add channels defaults
conda config --add channels r
conda config --add channels bioconda
conda update conda
安装配置软件：
-sratoolkit 教程上用的Linux方法，mac会报错，建议使用以下方法。 方法一（别人可以，我的不行）： conda install -c bioconda sra-tools=2.8.1
方法二（亲测可行）： ruby -e “$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)”
brew install wget
brew tap homebrew/science
brew install sratoolkit
conda install fastqc
conda install hisat2
conda install samtools
conda install -c bioconda htseq
conda install r
conda install rstudio
参考文献http://www./thread-1796-1-1.html http://www./thread-1800-1-1.html
读文献下数据AKAP95 regulates splicing through scaffolding RNAs and RNA processing factors. Nat Commun 2016 Nov 8;7:13347. PMID: 27824034 找数据地址  去GSE网站，搜GSE81916   下载代码
# 在终端运行
#用axel下载会快一些
brew install axel
for i in {56..62}
do
# 也可以用wget 下载
axel ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP075/SRP075747/SRR35899$i/SRR35899$i.sra
done
参考文献：1. http://www./thread-1908-1-1.html
了解fastqc需要用安装好的sratoolkit把sra文件转换为fastq格式的测序文件，并且用fastqc软件测试测序文件的质量！ 作业，理解测序reads，GC含量，质量值，接头，index，fastqc的全部报告，搜索中文教程，并发在论坛上面。 来源于生信技能树：http://www./forum.php?mod=viewthread&tid=1750#lastpost 实验步骤1. 将 sra 数据转化成 fastq 格式先建立一个SRR_fastqc.sh 文件，写入代码 #!/usr/bin/env bash
for i in {56..62}
do
fastq-dump --gzip --split-3 -O /Users/chengkai/Desktop/zhuanlu_files -A SRR35899${i}.sra
done
2. 在终端运行# 这一步很慢，我跑了2小时，泡杯咖啡，欣赏一部电影吧
$ bash SRR_fastqc.sh
3. fastqc 检测测序文件质量创建一个文件夹 mkdir fastqc/
创建一个fastqc.sh脚本，写入如下代码 #!/usr/bin/env bash
fastqc -o ./fastqc/ -t 8 SRR3589956_1.fastq.gz SRR3589956_2.fastq.gz
fastqc -o ./fastqc/ -t 8 SRR3589957_1.fastq.gz SRR3589957_2.fastq.gz
fastqc -o ./fastqc/ -t 8 SRR3589958_1.fastq.gz SRR3589958_2.fastq.gz
fastqc -o ./fastqc/ -t 8 SRR3589959_1.fastq.gz SRR3589959_2.fastq.gz
fastqc -o ./fastqc/ -t 8 SRR3589960_1.fastq.gz SRR3589960_2.fastq.gz
fastqc -o ./fastqc/ -t 8 SRR3589961_1.fastq.gz SRR3589961_2.fastq.gz
fastqc -o ./fastqc/ -t 8 SRR3589962_1.fastq.gz SRR3589962_2.fastq.gz
bash fastqc.sh
4. 质量解读html 格式用浏览器打开 基本信息每个read各位置碱基的测序质量横轴碱基的位置（1-51），纵轴是质量分数，20表示1%的错误率，30表示0.1% 红色线代表中位数，蓝色代表平均数，黄色是25%-75%区间，触须是10%-90%区间 Warning 报警 如果任何碱基质量低于10,或者是任何中位数低于25 Failure 报错 如果任何碱基质量低于5,或者是任何中位数低于20
偏离度- 横轴碱基的位置（1-51） - 纵轴是tail的Index编号 - 检查reads中每一个碱基位置在不同的测序小孔之间的偏离度，蓝色代表偏离度小，质量好，越红代表偏离度越大，质量越差。 - 这个图主要是为了防止，在测序过程中，某些tail受到不可控因素的影响而出现测序质量偏低 reads质量的分布
GC 含量统计- 横轴碱基的位置（1-51） - 纵轴是碱基含量百分比 - 图中四条线代表A T C G在每个位置平均含量 - 当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。 - 本结果前10个位置，每种碱基频率有明显的差别，说明有污染。 - 当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"
序列平均GC含量分布图- 横轴是百分比； 纵轴是每条序列GC含量对应的数量 - 蓝色的线是程序根据经验分布给出的理论值，红色是真实值，两个应该比较接近才比较好 - 当红色的线出现双峰，基本肯定是混入了其他物种的DNA序列 - 偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL" 各位置N的reads比率- 当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”，统计N的比率 - 正常情况下，N值非常小 - 当任意位置的N的比例超过5%，报"WARN"；当任意位置的N的比例超过20%，报"FAIL"
reads 长度分布
统计不同拷贝数的reads的频率- 横坐标是duplication的次数，纵坐标是duplicated reads的数目,以unique reads的总数作为100% - 测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在 - 当非unique的reads占总数的比例大于20%时，报"WARN"；当非unique的reads占总数的比例大于50%时，报"FAIL"
接头含量- 此图衡量的是序列中两端adapter的情况 - 如果在当时fastqc分析的时候-a选项没有内容，则默认使用图例中的四种通用adapter序列进行统计 - 本例中adapter都已经去除，如果有adapter序列没有去除干净的情况，在后续分析的时候需要先使用cutadapt软件进行去接头
重复短序列这个图统计的是，在序列中某些特征的短序列重复出现的次数 我们可以看到1-8bp的时候图例中的几种短序列都出现了非常多的次数，一般来说，出现这种情况，要么是adapter没有去除干净，而又没有使用-a参数；要么就是序列本身可能重复度比较高，如建库PCR的时候出现了bias 对于这种情况，我的办法是可以cut掉前面的一些长度，可以试着cut 1bp
参考文献1. http://fbb84b26./share/s/3XK4IC0cm4CL22pU-r1HPcQQ2irG2836uQYm2iZAyh1Zwf3_ （青山屋主） 2. www./thread-2034-1-1.html (laofuzi) 3. http://www.jianshu.com/p/14fd4de54402 （lxmic） 4. https://zhuanlan.zhihu.com/p/20731723 (孟浩巍)
了解参考基因和注释在UCSC下载hg19参考基因组，从gencode数据库下载基因注释文件，并且用IGV去查看你感兴趣的基因的结构，比如TP53,KRAS,EGFR等等。 作业，截图几个基因的IGV可视化结构！还可以下载ENSEMBL，NCBI的gtf，也导入IGV看看，截图基因结构。了解IGV常识。 转录组文章分析思路研究结果第一部分：转录组测序的实验设计中一般都会有一部分生理生化的实验结果，可能是电镜图片，可能是生理指标的变化，可能是有效物质的含量。 研究结果第二部分：根据生理生化指标选择了几个样本进行转录组测序？测了4G,6G数据量? 测序质量是否良好？检测到3000个基因表达？样本特异表达的基因有多少？ 研究结果第三部分：差异分析组合间有多少差异基因？不同差异组合共有和特有的基因数量有多少？不同样本表达模式的差异怎样？ 研究结果第四部分：差异分析组合间差异基因主要是那些基因家族？转录因子？KEGG通路？GO分类？这些研究结果需要去解释这生理生化指标的变化的原因
下载参考基因组
# 下载文件
axel http://hgdownload.soe./goldenPath/hg19/bigZips/chromFa.tar.gz
# 解压文件
tar -zxvf chromFa.tar.gz
#解压后可以发现，参考序列是按照染色体号分开列出的，我们还需要把所有的序列写入到一个文件中。
cat *.fa > hg19.fa
#最后删除其他无用的文件
rm chr*.fa
下载注释文件
# GTF格式主要是用来描述基因的注释
axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/GRCh37_mapping/gencode.v27lift37.annotation.gtf.gz
# GFF文件是一种用来描述基因组特征的文件，现在我们所使用的大部分都是第三版）（GFF3）
axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/GRCh37_mapping/gencode.v27lift37.annotation.gff3.gz
# 解压并删除原来的文件
gzip -d gencode.v27lift37.annotation.gtf.gz
gzip -d gencode.v27lift37.annotation.gff3.gz
IGV软件的下载和安装
IGV软件的使用窗口
主窗口布局： tool bar（工具栏），menu bar（菜单栏），pop-up menus（弹出式菜单） 染色体上的红色盒子表示显示这部分染色体，显示完整染色体是红框会消失 尺度显示了染色体的可见部分，刻度线显示了染色体的位置，跨度列表显示了当前显示的碱基的数量 IGV在水平行显示的数据称为tracks。通常,每个tracks代表一个样本或实验。这个例子展示了甲基化、基因表达、拷贝数，LOH和突变数据 IGV也显示某些特性,比如在tracks中的基因。默认情况下,IGV在一个面板显示数据，在另一个面板显示数据特性。拖放一个track名称，将一个track从一个面板移动到另一个地方 Track名称列在最左边面板。名字的易读性取决于 tracks的高度，例如，track越小，它的名字的可读性越小 属性名称被列在顶部的属性面板。彩色块代表属性值,每个独特的值被都有一个独特的颜色。鼠标放在一个颜色块的附近来查看其属性值 8.
导入参考基因组及注释信息，查看感兴趣基因的结构
-现在可以导入sort 后的gtf 文件了。
参考文献http://www.jianshu.com/p/48b5a0972301 (GTF/GFF文件的差异及其相互转换) http://www.jianshu.com/p/3e545b9a3c68 （hoptop） http://fbb84b26./share/s/3XK4IC0cm4CL22pU-r1HPcQQ1lZQRc2nKQhn2SthRW24I8CZ （greenhillman）
序列比对下载 indexaxel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar -zxvf hg19.tar.gz #解压
正式比对命令行目录：Desktop/zhuanlu_files hg19目录：Desktop/zhuanlu_files／hg19 RNA-Seq Data: Desktop/zhuanlu_files/RNA-Seq 注释： -t 记录时间 -x hg19(index)文件路径 -1 -2 测序的两个fastq文件 -S 比对结果输出路径
mkdir -p RNA-Seq/aligned
for i in `seq 56 58`
do
hisat2 -t -x hg19/genome -1 RNA-Seq/SRR35899${i}.sra_1.fastq.gz -2 RNA-Seq/SRR35899${i}.sra_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam
done
SAMtoolsSAMtools 详解： http://www.jianshu.com/p/15f3499a6469 SAM 数据格式是目前高通量测序中存放比对数据的标准格式，但是SAM 文件都很大，非常占空间，所以需要转到bam文件，而用的就是SAMtools软件。 for i in `seq 56 58`
do
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
samtools index SRR35899${i}_sorted.bam
done
IGV 可视化
reads 计数实现这个功能的软件也很多，还是烦请大家先自己搜索几个教程，入门请统一用htseq-count，对每个样本都会输出一个表达量文件。 需要用脚本合并所有的样本为表达矩阵。参考：生信编程直播第四题：多个同样的行列式文件合并起来。 对这个表达矩阵可以自己简单在excel或者R里面摸索，求平均值，方差。看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等。
分析前解读原文中人类293cell的数据只有3个（SRR3589956-58）, 缺少一个对照重复。所以后续我只分析小鼠的四个样品（SRR3589959-62） 理论基础因表达定量有以下三个水平： 基因水平 转录水平 外显子水平
标准化 reads计数后得到的基因定量结果（count matrix），在进行不同维度的比较时需要进行不同的处理 比较同一个样本(within-sample)不同基因之间的表达情况---主要考虑转录本长度的影响 比较不同样本（across-sample）同一个基因的表达情况---主要考虑测序深度的影响 标准化的算法：RPKM（SE）, FPKM（PE），TPM, TMM，RSEM等等。 标准化之后才能进行差异表达分析
1. 下载小鼠index 文件和基因组注释文件axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz
2. 序列与index比对for i in `seq 59 62`
do
hisat2 -t -x mouse/mm10/genome -1 RNA-Seq/SRR35899${i}.sra_1.fastq.gz -2 RNA-Seq/SRR35899${i}.sra_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam
done
3. 用reads 排序bam 文件for i in `seq 59 62`
do
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam
samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam
done
4. HTSeq 进行基因水平的定量HTSeq 详解基本用法
htseq-count [options] <alignment_file> <gff_file>
: 比对到基因组后得到的SAM文件 （SAMtools包含一些perl 脚本可以将大多数的比对文件转换成SAM格式 ），注意基因组mapping时一定要用支持剪接的比对软件（splicing-aware aligner）进行比对软件如TopHat. HTSeq-count 需要用到SAM格式中的 CIGAR 区域的信息。 想要通过标准输入来传入 基因组mapping得到SAM文件，用 – 替换 即可 如果你是双端测序，必须要对SAM进行排序（单端可不必排序，但这里我也推荐对单端测序结果排序已减少内存消耗并提高软件效率），对read name或 染色体位置 进行排序皆可（这里我推荐按read name排序，因为通过位置排序我遇到过错误）。具体需要通过 -r 参数指定，所以排序请详见参数 -r : 包含单位信息的gff/GTF文件（gff文件格式），大多数情况下就是指注释文件; 由于GTF文件其实就是gff文件格式的变形，在这里同样可以传入GTF格式文件。
模型参数-f <format>, –format=<format> ：
-r <order>, –order=<order>
samtools sort -n accepted_hits_unique.bam -o accepted_hits_unique_name_sorted.bam # -n 按read name 排序 ，如果不指定则按染色体位置排序
-s <yes/no/reverse>, –stranded=<yes/no/reverse>
For stranded=no, a read is considered overlapping with a feature regardless of whether it is mapped to the same or the opposite strand as the feature. For stranded=yesand single-end reads, the read has to be mapped to the same strand as the feature. For paired-end reads, the first read has to be on the same strand and the second read on the opposite strand. For stranded=reverse, these rules are reversed. -a <minaqual>, –a=<minaqual>
-t <feature type>, --type<feature type>
-i <id attribute>, --idattr=<id attribute>
-m <mode>, --mode=<mode>
-o <samout>, --samout=<samout>
-q, --quiet
-h, --help
5. reads 计数# 安装HTSeq
conda install HTSeq
# 语法
Usage: htseq-count [options] <sam_file> <gff_file>
mkdir -p RNA-Seq/matrix/
htseq-count RNA-Seq/aligned/SRR3589959_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589959.count
htseq-count RNA-Seq/aligned/SRR3589960_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589960.count
htseq-count RNA-Seq/aligned/SRR3589961_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589961.count
htseq-count RNA-Seq/aligned/SRR3589962_count.sam mouse/gencode.vM10.annotation.gtf > RNA-Seq/matrix/SRR3589962.count
6. 合并表达矩阵options(stringsAsFactors = FALSE)
# 首先将四个文件分别赋值：control1，control2，rep1，rep2
control1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))
control2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2"))
rep1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951"))
rep2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))
# 将四个矩阵按照gene_id进行合并，并赋值给raw_count
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(rep1,rep2, by="gene_id"))
# 需要将合并的raw_count进行过滤处理，里面有5行需要删除的行，在我们的小鼠的表达矩阵里面，是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter
raw_count_filt <- raw_count[-48823:-48825, ]
raw_count_filter <- raw_count_filt[-1:-2, ]
# 因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。
# 第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filter$gene_id)
# 将ENSEMBL重新添加到raw_count_filter矩阵
row.names(raw_count_filter) <- ENSEMBL
# 看一些基因的表达情况，在UniProt数据库找到AKAP95的id，并从矩阵中找到访问，并赋值给AKAP95变量
AKAP95 <- raw_count_filter[rownames(raw_count_filter)=="ENSMUSG00000024045",]
# 查看AKAP95
AKAP95
7. 简单分析8. 参考文献http://www.jianshu.com/p/e9742bbf83b9 （hoptop 大神） http://www.jianshu.com/p/24cf44b610a7 （lxmic 大神） http://fbb84b26./share/s/3XK4IC0cm4CL22pU-r1HPcQQ3ouv263iAk012b-9tU2UXRgt（青山屋主大神）
差异基因分析1. 读取表达矩阵# 首先将四个文件分别赋值：control1，control2，rep1，rep2
control1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))
control2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2"))
rep1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951"))
rep2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))
# 将四个矩阵按照gene_id进行合并，并赋值给raw_count
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(rep1,rep2, by="gene_id"))
# 需要将合并的raw_count进行过滤处理，里面有5行需要删除的行，在我们的小鼠的表达矩阵里面，是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter
raw_count_filt <- raw_count[-48823:-48825, ]
raw_count_filter <- raw_count_filt[-1:-2, ]
# 因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。
# 第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filter$gene_id)
row.names(raw_count_filter) <- ENSEMBL
raw_count_filter <- raw_count_filter[ ,-1]
2. DESeq2# 下载DESeq2
source("https:///biocLite.R")
biocLite("GenomeInfoDb")
DEseq2要求输入数据是由整数组成的矩阵。 DESeq2要求矩阵是没有标准化的。
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #~在R里面用于构建公式对象，~左边为因变量，右边为自变量。
dds <- DESeq(dds) #标准化
res <- results(dds, contrast=c("condition","treated","control")) #差异分析结果
3. 构建dds矩阵表达矩阵：即上述代码中的countData，就是我们前面通过read count计算后并融合生成的矩阵，行为各个基因，列为各个样品，中间为计算reads或者fragment得到的整数。 样品信息矩阵：即上述代码中的colData，它的类型是一个dataframe（数据框），第一列是样品名称，第二列是样品的处理情况（对照还是处理等），即condition，condition的类型是一个factor。 差异比较矩阵： 即上述代码中的design。 差异比较矩阵就是告诉差异分析函数是要从要分析哪些变量间的差异，简单说就是说明哪些是对照哪些是处理。
# 这一步很关键，要明白condition这里是因子，不是样本名称；小鼠数据有对照组和处理组，各两个重复
> condition <- factor(c(rep("control",2),rep("akap95",2)), levels = c("control","akap95"))
# 获取count数据
> countData <- raw_count_filter[,1:4]
> colData <- data.frame(row.names=colnames(raw_count_filter), condition)
> dds <- DESeqDataSetFromMatrix(countData, colData, design= ~ condition)
# 查看一下dds的内容
> head(dds)
4. DESeq标准化dds# normalize 数据
> dds2 <- DESeq(dds)
# 查看结果的名称，本次实验中是 "Intercept"，"condition_akap95_vs_control"
> resultsNames(dds2)
# 将结果用results()函数来获取，赋值给res变量
res <- results(dds2)
# summary一下，看一下结果的概要信息
summary(res)
5. 提取差异分析结果# 获取padj（p值经过多重校验校正后的值）小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
> table(res$padj<0.05)
> res <- res[order(res$padj),]
> diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
> diff_gene_deseq2 <- row.names(diff_gene_deseq2)
> resdata <- merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)
# 得到csv格式的差异表达分析结果
> write.csv(resdata,file= "control_vs_akap95.cvs",row.names = F)
6. 参考文献https://zhuanlan.zhihu.com/p/30350531（青山屋主） http://www.jianshu.com/p/3bfb21d24b74 （lxmic） http://www./packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html （DESeq2 官网）
|
