RNA-seq基本分析流程
1.从NCBI下载原始数据
多线程同时下载
Copy dataSRRID=(SRR8089897
SRR8089896
SRR8089895
)
for id in ${dataSRRID[@]};
do
j=`echo ${id:0:6}`
wget -c ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/${j}/${id}/${id}.sra &
done 检查数据大小
Copy ### 使用curl命令访问NCNI网页
for i in `ls .|grep sra`; do j=`echo ${i:0:10}`; curl https://www.ncbi.nlm.nih.gov/sra/${j}|grep -E "<tbody><tr>.*</tr></tbody>" >>1 ; done
### 正则表达式将数据大小给匹配出来，之后网页内容可能会有所改变，正则表达式可能不固定
sed 's/.*<tbody><tr>\(.*\)<\/tr><\/tbody>/\1/g' 1|sed 's/.*align=\"right\">\(.*\)<\/td><td>.*/\1/g' >2
### 与本地数据进行比较
ls -lh|grep sra |paste - 2 |less 将sra文件转换成fastq文件
Copy dataSRRID=(SRR8089897
SRR8089896
SRR8089895
)
for i in ${dataSRRID[@]};
do
fastq-dump --split-3 ./rawdata/${i}.sra -O ./fast_dump/ &
done 2.数据过滤 Trimmomatic
代码如下
Copy java -jar ~/software/Trimmomatic-0.39/trimmomatic-0.39.jar
PE
-threads 10
R1.fastq.gz
R2.fastq.gz
-baseout 输出文件目录/文件前称
ILLUMINACLIP:/public/home/zpliu/software/Trimmomatic-0.39/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:10 TRAILING:10 SLIDINGWINDOW:4:20 MINLEN:50 R1.fastq.gz/R2.fastq.gz 对应双端测序数据
baseout 后面接输出文件路径，以及输出文件的前缀
LEADING 去除头部read质量低于10的序列
SLIDINGWINDOW 按照4个碱基长度进行滑动，去除平均质量低于20的read
写一个bash循环进行批量过滤原始数据
Copy for i in ${A2list[@]};do
{
A2b=${i/R1/R2}
A2d=${i/R1/}
java -jar ~/software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 10 ${c}/A2/${i} ${c}/A2/${A2b} -baseout ${c}/A2_Trimmomatic/${A2d} ILLUMINACLIP:/public/home/zpliu/software/Trimmomatic-0.39/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:10 TRAILING:10 SLIDINGWINDOW:4:20 MINLEN:50
} &
done &作用是将for循环放入后台，进行并行计算，基本上100多G的数据 30分钟跑完
3.将数据比对到参考基因组
构建参考基因组索引'
接收基因组fasta序列文件，和索引生成路径及索引前称
Copy echo "Ghirsutum_genome_HAU_v1 begin build"
hisat2-build /public/home/zpliu/genome_data/genome_Ghitsutum_NAU/genome.Ghir.NAU.fa ./genome_Ghitsutum_NAU/Ghitsutum 进行比对
Copy Gh_indexfile='/public/home/zpliu/Hisat2Index/Ghirsutum_genome_HAU_v1.1/Ghirsutum_genome_HAU_v1'
Gh_hisatout="./hisat2out/"
all_fastq=`ls ./Trimmomatic|grep "_1P"`
for i in ${all_fastq[@]};
do
j=`echo ${i}|sed 's/_1P/_2P/g'`
k=`echo ${i}|sed 's/_1P//g'`
hisat2 -x ${Gh_indexfile} -1 ./Trimmomatic/${i} -2 ./Trimmomatic/${j} -p 10 --known-splicesite-infile /public/home/zpliu/genome_data/Ghirsutum_genome_HAU_v1.1/hista_splice.txt -S ./hisat2out/${k}.sam && samtools view -S ./hisat2out/${k}.sam -@ 10 -b -o ./hisat2out/${k}.bam && samtools sort -@ 10 ./hisat2out/${k}.bam -O bam -o ./hisat2out/${k}_sort.bam
done 4.计算基因表达量
将比对好的sam文件按照染色体位置进行排序后，使用stringtie比对
Copy Gh_gff="/public/home/zpliu/genome_data/Ghirsutum_genome_HAU_v1.1/Ghirsutum_gene_model.gff3"
for i in `ls ./hisat2out`;
do
j=`echo ${i}|sed 's/_.*//g'`
stringtie ./hisat2out/${i} -G ${Gh_gff} -e -p 10 -A ./stringtie/${j}
done 合并多个组织的表达量
Copy # 构造字典文件
sed 's/^/[/g' 111 |sed 's/$/]=/g'|paste - 22 |sed 's/\t//g' >end
/*
[SRR8090044]=TM1_10DPA_Fiber_1
[SRR8090041]=TM1_10DPA_Fiber_2
[SRR8090042]=TM1_10DPA_Fiber_3
[SRR8090046]=TM1_15DPA_Fiber_1
[SRR8090049]=TM1_15DPA_Fiber_2
[SRR8090050]=TM1_15DPA_Fiber_3
[SRR8090004]=TM1_20DPA_Fiber_1
[SRR8090007]=TM1_20DPA_Fiber_2
[SRR8090006]=TM1_20DPA_Fiber_3
[SRR8090010]=TM1_25DPA_Fiber_1
*/
declare -A sample
sample=(end文件中的内容)
# 对每个样本的数据盖头换面
sort SRR8090089|cut -f 8|sed '1d'|awk 'BEGIN{print "'$a'"}{print $0}'|less
#SRR数组可以按照一定顺序输出
SRR=(SRR8090087 SRR8090088)
# 使用字典批量,这里由于每个生成文件中基因数目不一致，先提取共有的基因存进geneid里，并行运算
for i in ${SRR[@]};
do {
cat geneid|xargs -I {} grep {} $i |cut -f 8|awk 'BEGIN{print "'${sample[$i]}'"}{print $0}' >${i}_tmp; } &
done 手动计算指定区域的表达量
链特异性建库
使用samtools 计算目标区域的表达量
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
需要事先利用利用上面讲的sort和建index命令执行完后，用下述命令。
Usage: samtools tview [ref.fasta]
Copy samtools tview -p Gbar_A07:41137452 align_sorted.bam genome.fa 出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。 按下 g ，则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第 10号scaffold的第1000个碱基位点处。 使用H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。 使用空格建向左快速移动（和 L 类似），使用Backspace键向左快速移动（和 H 类似）。 Ctrl+H 向左移动1kb碱基距离； Ctrl+L 向右移动1kb碱基距离 可以用颜色标注比对质量，碱基质量，核苷酸等。30～40的碱基质量或比对质量使用白色表示； 20～30黄色；10～20绿色；0～10蓝色。 使用点号'.'切换显示碱基和点号；使用r切换显示read name等 还有很多其它的使用说明，具体按 ？ 键来查看。
Copy 结果说明：
“.” 比对到正链;
“，” 表示比对到负链;
“<”或“>” 表示reference skip RNA-seq当中内含子剪切;
"ATCGN" 表示正向mismatch;
"atcgn" 表示反向mismatch;
‘+[0-9]+[ACGTNacgtn]+’ insertion;
‘-[0-9]+[ACGTNacgtn]+’ 表示deletion;
“^”标记reads起始;
“$”标记reads segment结尾; 计算某个区域paired read数目
Copy ##提取对应区域的bam文件
samtools view -O BAM align_sorted.bam Gbar_A07:41137452-41138228 >Gbar_A07G016470.bam
##按照read名字排序后，转成bed文件
samtools sort -n Gbar_A07G016470.bam -O BAM -o Gbar_A07G016470_sortReadName.bam
##提取对应的bed文件
bamToBed Gbar_A07G016470_sortReadName.bam -bedpe -mate1 >11.bed
#如果不存在成对比对的read，则会跳过当前read
其中最后两列表示insert 片段中read1 和read1比对到的链
+ RNA信息是非模板链（基因所在链的信息）
+ 反转录得到cDNA（模板链上的信息）
+ 经过dTUP处理后保留的是非模板链的信息
+ 再经过合成后read保留的就是模板链的信息
一般read1与基因的方向是相反的，而read2与基因方向是相同的
##区分基因组不同位置，正负链上比对到的read数目
+ 首先筛选比对质量大于20的read，且是成对的read
Copy ##获取read比对到基因组的正负链以及坐标
awk '$1!="."&&$8>=20{print $0}' 11.bed |awk '$10=="-"{print $1"\t"$2"\t"$6"\t"$10"\t"$7}$10=="+"{print $1"\t"$5"\t"$3"\t"$10"\t"$7}' >readcount.txt
##与基因坐标取交集
intersectBed -loj -a gene_1.bed -b readcount.txt |awk '$7=="+"{a[$2][1]+=1}$7=="-"{a[$2][2]+=1}$7=="."{a[$2][1]=0;a[$2][2]=0}END{for(i in a){print $1"\t"i"\t"a[i][1]"\t+\n"$1"\t"i"\t"a[i][2]"\t-"}}'|awk -F "\t" '$3==""{print $1"\t"$2"\t0\t"$4}$3!=""{print $0}'|sort -k2,2n >read_align_gene.txt 使用脚本计算某区域的read比对情况
脚本只支持链特异性建库的bam文件
Copy ##提取某个基因组区域比对到的read数目
python /public/home/zpliu/github/jupyter/ResearchProject/fiber-DRGE/script/sQTL/lncRNA_quality.py BamFile outFile sampleName Ghir_A07:42970060-42970318:+