##说明：此代码是运行处理depseq的程序

#文件位置 实验室大服务器
# 运行目录 /home/chaim/disk/depseq/

#文档结构 

#测序原始数据
# R28_cleaned_R1.fastq.gz  R29_cleaned_R1.fastq.gz  r28_cleaned_R1.fastq.gz  r29_cleaned_R1.fastq.gz
# R28_cleaned_R2.fastq.gz  R29_cleaned_R2.fastq.gz  r28_cleaned_R2.fastq.gz  r29_cleaned_R2.fastq.gz

#zm437 /玉米基因组zm437版本序列

1.使用bowtie2比对

#1.1 建立基因组的索引文件
 bowtie2-build zm437 index
#1.2开始比对
bowtie2 -p 8 -x index -1 R28_cleaned_R1.fastq.gz -2 R28_cleaned_R2.fastq.gz -S R28.sam &
bowtie2 -p 8 -x index -1 R29_cleaned_R1.fastq.gz -2 R29_cleaned_R2.fastq.gz -S R29.sam 
#bowtie2 -p 8 -x index -1 r28_cleaned_R1.fastq.gz -2 r28_cleaned_R2.fastq.gz -S r28.sam &
#bowtie2 -p 8 -x index -1 r29_cleaned_R1.fastq.gz -2 r29_cleaned_R2.fastq.gz -S r29.sam

#1.2对比对的sam结果添加read group信息
bowtie2 -p 8 --rg-id R28 --rg "PL:ILLUMINA" --rg "SM:R28" -x index -1 R28_cleaned_R1.fastq.gz -2 R28_cleaned_R2.fastq.gz -S R28ad.sam 
bowtie2 -p 8 --rg-id R29 --rg "PL:ILLUMINA" --rg "SM:R28" -x index -1 R29_cleaned_R1.fastq.gz -2 R29_cleaned_R2.fastq.gz -S R29ad.sam 
bowtie2 -p 8 --rg-id r28 --rg "PL:ILLUMINA" --rg "SM:R28" -x index -1 r28_cleaned_R1.fastq.gz -2 r28_cleaned_R2.fastq.gz -S r28ad.sam 
bowtie2 -p 8 --rg-id r29 --rg "PL:ILLUMINA" --rg "SM:R28" -x index -1 r29_cleaned_R1.fastq.gz -2 r29_cleaned_R2.fastq.gz -S r29ad.sam

#1.3samtools (对sam文件进行排序，并转换Sam为bam)
samtools sort R28ad.sam > R28.bam &
samtools sort R29ad.sam > R29.bam &
samtools sort r28ad.sam > r28.bam &
samtools sort r29ad.sam > r29.bam 
#1.3.2 samtools 合并两个bam
samtools cat R28.bam R29.bam -o R.bam
samtools cat r28.bam r29.bam -o r.bam
#1.4 samtools 建立索引
samtools index R.bam
samtools index r.bam

2.macs2 构建表达峰图

参考地址1
参考地址2

#macs2 callpeak -t R29.bam -n sample --shift -100 --extsize 200 --nomodel -B --SPMR -g hs --outdir R29_out 2> sample.macs2.log
macs2 callpeak -t R.bam -c r.bam -g 2.1e109 -B -f BAMPE -n R -q 0.00001

参数讲解：
-f BAMPE 双端测序
-g hs 基因组实际大小，此处是人类的hg,如果没有对应的物种信息，需要自己设定大小。
3.snpEFF
注释突变位点信息

java -Xmx128g -jar /home/chaim/bsa/snpEff/snpEff/snpEff.jar -c  /home/chaim/bsa/snpEff/snpEff/snpEff.config zm437 ${name_array[$n]}"_peaks.xls" > ${name_array[$n]}".snp.eff.vcf" -stats ${name_array[$n]}".html"

4.bedtools
主要是bedtools转换数据格式

先使用grep 删除vcf文件头部的内容，只保留从染色体1开始的数据信息。前面的头部文件一定要删除，否则后续bedtools无法完成转换

#查看AY.snp.eff.vcf有多少行以#开头
      grep -n '#' AY.snp.eff.vcf  |wc -l  
      grep -n '#' GY.snp.eff.vcf |wc -l
#删除前面多余的行，此处是31行，删除前面的31行
      sed '1,31d' AY.snp.eff.vcf  >AY.eff.vcf
      sed '1,31d' GY.snp.eff.vcf >GY.eff.vcf

使用bedtools转换bed格式为fasta格式。

    bedtools getfasta -fi zm437.fa -bed AY.eff.vcf -fo AY.out.fa
    bedtools getfasta -fi zm437.fa -bed GY.eff.vcf -fo GY.out.fa

5.MEME搜寻motif