sra文件地址：ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/

下載測序數(shù)據(jù)

# 循環(huán)下載sra文件

1扇单，for ((i=677;i<=680;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP029/SRP029245/SRR957$i/SRR957$i.sra ;done

# 后臺運行,注意括號的位置

for ((i=508;i<=523;i++)) ;do(nohup wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR103/SRR1039$i/SRR1039$i.sra &) ;done

2，ncbi上獲得該項目的Accession List（如圖16個run代表16個數(shù)據(jù)）是一個txt文件焚虱，改名為id

image

cat id | while read id ;do prefetch $id ;done  ##  while read 一次讀取一行
ps -ef | grep prefetch  ## 查看正在運行的prefetch命令订框，ps -ef 相當于windows的任務管理器
ps -ef | grep prefetch | awk '{print $2}' | while read id ;do kill $id ;done ## awk提取出任務編號所在的第二列，循環(huán)殺除任務
ls | while read id ;do (nohup ../../biosoft/sratoolkit/bin/fastq-dump --gzip --split-3 -o ./ $id &) ;done  ## 循環(huán)把下載的所有sra文件都變?yōu)閒astq 

測序數(shù)據(jù)質(zhì)控

ls *.gz | while read id ;do fastqc $id ;done # 循環(huán)fastqc處理每個fastq文件
ls *.gz | xargs fastqc # 與上等效晨另，xargs將ls的輸出內(nèi)容作為參數(shù)傳遞給fastqc,一次傳遞一個
multiqc ./ # 把每個數(shù)據(jù)的fastqc質(zhì)控報告潭千，合并到一個報告里，方便查看
ls *1.fastq.gz >1 
ls *2.fastq.gz >2  
paste 1 2 > config # 把所有的1.fastq.gz和2.fastq.gz的文件名合并到config文件

#寫一個循環(huán)修剪的shell腳本借尿，命名為qc.sh刨晴，軟件為trim_galore

dir='/sas/supercloud-kong/liuhuijie/RNA-seq/clean'  # 輸出地址

cat $1 | while read id # $1 表示輸命令時的第一個參數(shù)

do

arr=$id

fq1=${arr[0]} # config文件里一行是兩個文件名，中間被空格隔開arr[0]表示第一個

fq2=${arr[1]}

nohup trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &

done

bash qc.sh config # 執(zhí)行命令

比對路翻，把所有的fastq數(shù)據(jù)狈癞，提取出前10000行，方便比對茂契。

1. ls ../*.gz | while read id ;do (zcat $id | head -10000 > $(basename $id)) ;done # 提取前1萬行蝶桶，重新輸出到$(basename $id)文件里，basename是只取文件名
2. ls ../*.gz | while read id ;do (zcat $id | head -10000 > $(basename $id '.gz')) ;done # 重輸出的文件名里去除了后綴 .gz

分別用軟件bowtie2掉冶、hisat2做比對得到sam文件真竖，批量轉(zhuǎn)換為bam文件

hisat2 -p 5 -x ../../index/hisat/hg38 -1 SRR1039508_1_val_1.fq -2 SRR1039508_2_val_2.fq -S tmp.hisat.sam
bowtie2 -p 5 -x ../../index/bowie/hg38 -1 SRR1039508_1_val_1.fq -2 SRR1039508_2_val_2.fq > tmp.bowtie.sam
ls *.sam | while read id ;do (samtools sort -O bam -@ 5 -o $(basename $id '.sam').bam $id) ;done  # 批量轉(zhuǎn)換為bam，把后綴.sam去掉換成bam
samtools view tmp.bowtie.bam | less -S  #查看bam文件厌小，可試試IGV可視化查看
ls *.bam | xargs -i samtools index {}  #批量為bam文件建索引恢共，xargs -i把bam文件作為參數(shù)傳遞給samtools，{}代替了bam文件
ls *bam | while read id ;do (samtools flagstat -@ 10 $id > $(basename $id '.bam').flagstat) ;done  # flagstat 統(tǒng)計
cat *flagstat | cut -d ' ' -f 1 | paste - - - - - - - - - - - - - # 合并所有的flagstat文件的第一列璧亚，每個flagstat文件都為13行讨韭，有13個-

把當前文件夾里的fastq數(shù)據(jù)，用hisat2和bowtie軟件癣蟋，按順序比對到hg38

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 5 -x ../index/hisat/hg38 -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S $id.hisat.sam);done

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 | bowtie2 -p 5 -x ../index/bowtie/hg38 -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz > $id.bowtie.sam);done

定量

nohup featureCounts -T 10 -p -t exon -g gene_id -a ../gtf/Homo_sapiens.GRCh38.92.chr.gtf.gz -o all.id.txt *bam 1>count.log 2>&1 &
#  featureCounts 進行定量透硝，統(tǒng)計比對在這個基因的坐標上的read數(shù)
multiqc all.id.txt.summary  #用multiqx對上步得到的結(jié)果進行可視化

表達矩陣探索

rm(list = ls())
options(stringsAsFactors = F)
a= read.table('all.id.txt',header = T)  #header
meta=a[,1:6]
exprSet=a[,7:ncol(a)]
a2=exprSet[,'SRR1039508.hisat.bam']  ##
library(pheatmap)
png('heatmap.png')
corrplot(cor(exprSet))
pheatmap(scale(cor(log2(exprSet+1))))
dev.off()
## 
library(airway)
data("airway")
exprSet=assay(airway)
a1=exprSet[,'SRR1039508']  ##
group_list=colData(airway)[,3]
##hclust, 層次聚類包
colnames(exprSet)=paste(group_list,1:ncol(exprSet),sep = '_')
##difine 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))  #par函數(shù)設(shè)置圖形邊距，mar參數(shù)設(shè)置邊距
png('hclust.png',res = 120)
plot(as.dendrogram(hc),nodePar=nodePar,horiz=TRUE)
dev.off()

a2=data.frame(id=meta[,1],a2=a2)
a1=data.frame(id=names(a1),a1)
library(stringr)  #stringr包處理字符串
a2$id=str_split(a2$id,'\\.',simplify = T)[,1]
tmp=merge(a1,a2,by='id')
png('tmp.png')
plot(tmp[,c(2,4)])  #從圖上看更直觀
dev.off()

DEseq2篩選差異表達基因

library(DESeq2)
library(edgeR)
library(limma)
library(airway)
data("airway")
exprSet=assay(airway)  #定量后的信息
group_list=colData(airway)[,3]  #提取出來了分組信息疯搅，也可以手動寫成c()
colData=data.frame(row.names = colnames(exprSet),
                   group=group_list)
dds=DESeqDataSetFromMatrix(countData = exprSet,
                           colData = colData,
                           design = ~group)  #獲取矩陣信息
dds=DESeq(dds)
res=results(dds,contrast = c('group','trt','untrt'))
summary(res)  #利用summary命令統(tǒng)計顯示一共多少個genes上調(diào)和下調(diào)
resOrdered=res[order(res$padj),]  #根據(jù)padj（p值經(jīng)過多重校驗校正后的值）排序
DEG=as.data.frame(resOrdered)
DEG=na.omit(DEG)
library(pheatmap)
choose_gene=head(rownames(DEG),100)  
choose_matrix=exprSet[choose_gene,]  #抽取差異表達顯著的前100個基因
choose_matrix=t(scale(t(choose_matrix)))  #用t函數(shù)轉(zhuǎn)置蹬铺，scale函數(shù)標準化
pheatmap(choose_matrix,filename='DEG_top100.png')

火山圖

logFC_cutoff=with(DEG,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))  #算log2FoldChange的閾值，with 提取數(shù)據(jù)框中的某些參數(shù)做運算秉撇，abs求絕對值，sd求標準差
DEG$change=as.factor(ifelse(DEG$pvalue<0.05 & abs(DEG$log2FoldChange)>logFC_cutoff,ifelse(DEG$log2FoldChange>logFC_cutoff,'UP','DOWN'),'NOT'))
#ifelse函數(shù)秋泄，大于logFC_cutoff的設(shè)為up琐馆，小于為down
this_title=paste0('Cutoff for logFC is ',round(logFC_cutoff,3), '\nThe number of up gene is ',nrow(DEG[DEG$change=='UP',]), '\nThe number of down gene is ',nrow(DEG[DEG$change=='DOWN',]))  #paste0函數(shù)，默認是sep=""
library(ggplot2)                     
g=ggplot(data = DEG,
         aes(x=log2FoldChange,y=-log10(pvalue),
             color=change))+
  geom_point(alpha=0.4,size=1.75)+
  theme_set(theme_set(theme_bw(base_size = 20)))+
  xlab('log2 fold change')+ylab('-log10 p-value')+
  ggtitle(this_title)+theme(plot.title = element_text(size = 15,hjust = 0.5))+
  scale_color_manual(values = c('blue','black','red'))  #corresponding to the levels(res$change)
ggsave(g,filename = 'volcano.png') 

png('dispersions.png',1000,1000,pointsize = 20)
plotDispEsts(dds,main='dispersion plot')
dev.off()

作圖查看原始定量后的數(shù)據(jù)和normalization后的數(shù)據(jù)的差異

rld=rlogTransformation(dds)  #DEseq2自帶的rlog算法對數(shù)據(jù)進行count矩陣轉(zhuǎn)換
exprMatrix_rlog=assay(rld)
png('DEseq_RAWvsNORM.png',height = 800,width = 800)
par(cex=0.7)  #par函數(shù)設(shè)定全局繪圖參數(shù)恒序，cex放大多少倍
n.sample=ncol(exprSet)
if(n.sample>40) par(cex=0.5)  #
cols=rainbow(n.sample*1.2)  #rainbow 漸變的彩虹色
par(mfrow=c(2,2))  #4個圖按行排序
boxplot(exprSet, col=cols,main='expression value',las=2)  #las為2瘦麸，標簽垂直坐標軸
boxplot(exprMatrix_rlog, col=cols, main='expression value',las=2)
hist(as.matrix(exprSet))
hist(exprMatrix_rlog)
dev.off()