转录组全新学习之总篇
2018-07-04 本文已影响44人
chaimol
- 数据已经存在服务器,基因组注释文件、基因组文件都已经存在。服务器环境软件都已经安装完成。但是没有root权限。
数据获取预处理
find ./510 -name '*zaoqian*'|xargs -i mv {} ./100 #批量查找510目录里所有包含zaoqian的文件,移动到100目录里。
- 实验流程和使用软件
- 质控 trimmomatic
trimmomatic功能讲解
精华篇 trimmomatic参数讲解 - 比对 HISAT2
精华HISAT2 - 装配 String tie
精华stringtie - 定量分析 Deseq2
- 转录本分析 Cufflink
- 从头组装 trinity
- 开始第一步——质控。
fastqc可以检测测序质量。
fastqc *.fastq.gz -t 4 #开启4个进程。
0.jpg
trimmomatic去除接头和低质量测序数据(-threads 12 开启了12个进程)
java -jar /home/guo/tool/Trimmomatic-0.38/trimmomatic-0.38.jar PE -threads 12 -phred33 CR-zaoqian-1_L8_I313.R1.clean.fastq.gz CR-zaoqian-1_L8_I313.R2.clean.fastq.gz CR-paired-1-R1.fastq.gz CR-unpaired-1-R1.fastq.gz CR-paired-1-R2.fastq.gz CR-unpaired-1-R2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 &
最后的&是在后台执行。
- 比对HISAT2
- 建立索引文件
#-p 8 8个进程,基因组文件Zea_mays.AGPv4.cds.all.fa /disks/···/maize437 指定输出的文件的名字为maize437,存放位置就是前面这一串路径。
hisat2 -build -p 8 Zea_mays.AGPv4.cds.all.fa /disks/backup/chaim/maize437
- 开始比对
hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 /disks/backup/chaim/cms/zaoqian/Cr-paired-1-R1.fastq.gz -2 /disks/backup/chaim/cms/zaoqian/Cr-paired-1-R2.fastq.gz -S Cr-1.map.sam --dta-cufflinks --novel-splicesite-outfile Cr-1.nsplice &
参数讲解
-x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel #基因组索引文件
-p 8 #使用8个进程
-1 paired1 #质控后左链的数据
-2 paired2 #质控后右链的数据
-S Cr-1.map.sam #输出文件
-dta-cufflinks--novel-splicesite-outfile Cr-1.nsplice #后期使用cufflinks做准备。
- stringtie官网文档(自带梯子)
for (( i=1;i<4;i++ ));
do
stringtie CR.map.bam -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o CR-${i}.gtf &
stringtie Nr.map.bam -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o Nr-${i}.gtf &
stringtie Cr.map.bam -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o Cr-${i}.gtf &
done
stringtie --merge -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -o merged.gtf CR-1.gtf CR-2.gtf CR-3.gtf Cr-1.gtf Cr-2.gtf Cr-3.gtf Nr-1.gtf Nr-2.gtf Nr-3.gtf &
for (( i = 0; i < 4; i++ ));
do
stringtie Cr-${i}.map.bam -G merged.gtf -p 8 -b ${A}${i}"_out" -e -o Cr-${i}ss.gtf &
stringtie Nr-${i}.map.bam -G merged.gtf -p 8 -b ${B}${i}"_out" -e -o Nr-${i}ss.gtf &
stringtie CR-${i}.map.bam -G merged.gtf -p 8 -b ${C}${i}"_out" -e -o CR-${i}ss.gtf &
done
参数讲解:
-G 参考基因组注释文件
-p 8个进程
-o 输出注释文件名
-b 输出其他文件的路径(会同时生成多个文件,一定要输出在不同的路径中,不然后面的输出会覆盖前面的输出。)
6.格式转换,数据整合。
使用的是python的脚本。下载地址
python2.7 /disks/backup/chaim/soft/prepDE.py -i gtf2
#注意此处的prepDE.py的python版本是python2.7,一定要指定版本号。否则会报错。
gtf2文件内容是对应的比对后文件名和文件的位置。
Cr1-st ./Cr1-st.gtf
Cr2-st ./Cr2-st.gtf
Cr3-st ./Cr3-st.gtf
CR1-st ./CR1-st.gtf
CR2-st ./CR2-st.gtf
CR3-st ./CR3-st.gtf
Nr3-st ./Nr3-st.gtf
Nr1-st ./Nr1-st.gtf
Nr2-st ./Nr2-st.gtf
输出文件为
转录本输出矩阵 transcript_count_matrix.csv
基因输出矩阵 gene_count_matrix.csv
[perpDE.py参数](http://ccb.jhu.edu/software/stringtie/index.shtml?
t=manual#deseq)
-i 输入文件
-g gene输出矩阵位置
-t 转录本输出矩阵位置
- Deseq2定量分析
安装在本地Ubuntu环境
使用conda安装
>conda install biocpnductor-deseq2
安装到服务器端
source("[http://bioconductor.org/biocLite.R](http://bioconductor.org/biocLite.R)")
biocLite("DESeq2", dependencies=T)
使用DESeq2
- Cufflinks
整个实验的大概脚本情况
#!/bin/bash
#批量循环,创造出文件名,并且复制文件。
#mkdir 102
A="CR"
B="Cr"
C="Nr"
tit[0]="_L8_I313." #CR-1
tit[1]="L8_I314." #CR-2
tit[2]="_L2_I301." #Cr-1
tit[3]="L2_I302." #Cr-2
tit[4]="_L2_I307." #Nr-1
tit[5]="L2_I308." #Nr-2
mid="-paired-"
end=".fastq.gz"
R1="-R1.fastq.gz"
R2="-R2.fastq.gz"
bam1="map.bam"
for (( i = 1; i < 4; i++ ));
do
k=$i
# echo ${A}${mid}${i}${R1}
# echo ${A}${mid}${i}${R2}
# echo ${B}${mid}${i}${R1}
# echo ${B}${mid}${i}${R2}
# echo ${C}${mid}${i}${R1}
# echo ${C}${mid}${i}${R2}
#第1步质控:trimmomatic(共9条,更改输入输出文件名即可)
#java -jar /home/guo/tool/Trimmomatic-0.38/trimmomatic-0.38.jar PE -threads 8 -phred33 Cr-zaoqian-1_L2_I301.R1.clean.fastq.gz Cr-zaoqian-1_L2_I301.R2.clean.fastq.gz Cr-paired-1-R1.fastq.gz Cr-unpaired-1-R1.fastq.gz Cr-paired-1-R2.fastq.gz Cr-unpaired-1-R2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 &
#第2步比对:hisat2(共三条)
#hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 ${A}${mid}${i}${R1} -2 ${A}${mid}${i}${R2} -S ${A}${i}".map.sam" --dta-cufflinks --novel-splicesite-outfile ${A}${i}".nsplice" &
#hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 ${B}${mid}${i}${R1} -2 ${B}${mid}${i}${R2} -S ${B}${i}".map.sam" --dta-cufflinks --novel-splicesite-outfile ${B}${i}".nsplice" &
#hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 ${C}${mid}${i}${R1} -2 ${C}${mid}${i}${R2} -S ${C}${i}".map.sam" --dta-cufflinks --novel-splicesite-outfile ${C}${i}".nsplice" &
#第3步:用samtool,格式转换,将sam转换为bam(共三条)
# samtools sort -@ 8 -o ${A}${i}".map.bam" ${A}${i}".map.sam" &
# samtools sort -@ 8 -o ${B}${i}".map.bam" ${B}${i}".map.sam" &
# samtools sort -@ 8 -o ${C}${i}".map.bam" ${C}${i}".map.sam" &
#第4步装配:用stringtie(共三轮)
#第一轮(9个分别比对到基因组)
#stringtie ${A}${i}".map.bam" -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o ${A}${i}".gtf" &
#stringtie ${B}${i}".map.bam" -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o ${B}${i}".gtf" &
#stringtie ${C}${i}".map.bam" -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o ${C}${i}".gtf" &
#第二轮(整合9个的结果成一个)
#stringtie --merge -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o merged.gtf CR1.gtf CR2.gtf CR3.gtf Cr1.gtf Cr2.gtf Cr3.gtf Nr1.gtf Nr2.gtf Nr3.gtf &
#第三轮(以第二轮的结果作为参考序列,9个分别比对)
# mkdir ${A}${i}"_out"
# mkdir ${B}${i}"_out"
# mkdir ${C}${i}"_out"
# stringtie ${A}${i}".map.bam" -G merged.gtf -p 8 -b ${A}${i}"_out" -e -o ${A}${i}"-st.gtf" &
# stringtie ${B}${i}".map.bam" -G merged.gtf -p 8 -b ${B}${i}"_out" -e -o ${B}${i}"-st.gtf" &
# stringtie ${C}${i}".map.bam" -G merged.gtf -p 8 -b ${C}${i}"_out" -e -o ${C}${i}"-st.gtf" &
第5步 生成CSV文件
#python路径 python2.7 /disks/backup/chaim/soft/prepDE.py -i
第6步 deseq2进行定量分析
done