精华文章收藏转录组分析走进转录组

转录组全新学习之总篇

2018-07-04  本文已影响44人  chaimol
  1. 数据已经存在服务器,基因组注释文件、基因组文件都已经存在。服务器环境软件都已经安装完成。但是没有root权限。
    数据获取预处理
find ./510 -name '*zaoqian*'|xargs -i mv {} ./100    #批量查找510目录里所有包含zaoqian的文件,移动到100目录里。

参考
1
2

  1. 实验流程和使用软件
  1. 开始第一步——质控。
    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 &

最后的&是在后台执行。

  1. 比对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做准备。

  1. 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 转录本输出矩阵位置

  1. Deseq2定量分析
    安装在本地Ubuntu环境
    使用conda安装
    >conda install biocpnductor-deseq2

安装到服务器端
source("[http://bioconductor.org/biocLite.R](http://bioconductor.org/biocLite.R)")
biocLite("DESeq2", dependencies=T)

使用DESeq2


  1. 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
上一篇下一篇

猜你喜欢

热点阅读