参考： http://192.168.47.160:8080/wiki/GB/doku.php?id=products_tr:pipelines:海外项目执行文档-wgbs

/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Pipeline/Methylation_Pipeline_V7.3.2.qsubname.py

<color #ed1c24>注意</color>：物种为人时，基因组推荐选择primary.fa（3G），如有关联分析或客户要求选择toplevel.fa（60G），提醒运营周期延长，并修改mapping内存投递为80G。repeat可直接使用临近版本已准备文件，cxx_prepare耗时较长可先进行QC。raw data进行拆分，周期紧张时DMR分析可进一步拆分。

ref=/TJPROJ13/GB_TR/reference_data/Animal/Mus_musculus/Mus_musculus_GRCm39/WGBS
python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Pipeline/Methylation_Pipeline_V7.3.2.qsubname.py \
        --project XXXXXXXXXXXXXX-Z01-J001_Mus_musculus,XXXXXX（委托）合同,XXXXXXXXXXXXXX \
        --sampleconfig samplelist.config \
        --HEADCROP 10 \
        --library_type LiBS \
        --species mmu \
        --direction yes \
        --datasize 100 \
        --group A1:A2:A3,B1:B2:B3 \
        --groupname A,B \
        --compare 2:1 \
        --autogenome /TJPROJ13/GB_TR/PJ_AI/AI_genome/plant/ensembl_42_prunus_persica_prunus_persica_ncbiv2_toplevel 如果填写了autogenome，则genome、gtf、goann、genenamefile、chrlist均不用填写
        --genome $ref \
        --gtf $ref/genome.gtf \
        --goann $ref/go.txt \
        --genenamefile $ref/genenamefile.txt \
        --chrlist $ref/chrlist.txt \
        --features CGI,CGI_shore,promoter,utr5,exon,intron,utr3,repeat \
        --name 生信名字 \
        --yunyingname 运营名字 \
        --dml no \
        --Profit gn \
        --context ALL \
        --clean N \
        --Data_Volume 90 \
        --queue gball1s.q,gbtr2s.q,gball1m.q,gbtr2m.q

python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Pipeline/Methylation_Pipeline_V7.3.2.qsubname.py \
        --project XXXXXXXXXXXXXX-Z01-J001_Mus_musculus,XXXXXX（委托）合同,XXXXXXXXXXXXXX \
        --sampleconfig samplelist.config \
        --HEADCROP 0 \
        --library_type RRBS \
        --species mmu \
        --direction yes \
        --datasize 100 \
        --group A1:A2:A3,B1:B2:B3 \
        --groupname A,B \
        --compare 2:1 \
        --autogenome /TJPROJ13/GB_TR/PJ_AI/AI_genome/plant/ensembl_42_prunus_persica_prunus_persica_ncbiv2_toplevel 如果填写了autogenome，则genome、gtf、goann、genenamefile、chrlist均不用填写
        --genome $ref \
        --gtf $ref/genome.gtf \
        --goann $ref/go.txt \
        --genenamefile $ref/genenamefile.txt \
        --chrlist $ref/chrlist.txt \
        --features promoter,utr5,exon,intron,utr3,CGI,CGI_shore,repeat \
        --name 生信名字 \
        --yunyingname 运营名字 \
        --dml no \
        --dmr_smoothing FALSE \  #有生物学重复
        --Profit gn \
        --context CG \
        --clean N \
        --Data_Volume 10 \
        --queue gball1s.q,gbtr2s.q,gball1m.q,gbtr2m.q

注意：RRBS比LiBS/WGBS多一行：–dmr_smoothing FALSE；–context填CG。

• 表中加粗参数均为必填参数
• 做差异分析及相应功能富集可选参数：–group，–groupname,–compare,–species,–genenamefile
  
• 表2为可选参数
• –cut，截取数据时填写
• –scoremin,调整mapping参数时填写

hw_LiBS:
--HEADCROP 10, --library_type LiBS, --dml no, --Profit hw 

hw_WGBS:
--HEADCROP 0, --library_type WGBS, --dml no, --Profit hw 

gn_LiBS: 
--HEADCROP 10, --library_type LiBS, --dml no, --Profit gn 

hw_RRBS 
--HEADCROP 0, --library_type RRBS, --dml no, --context CG, --Profit hw 

gn_RRBS 
--HEADCROP 0, --library_type RRBS, --dml no, --context CG, --Profit gn

作用：根据lims里的项目信息（final_xinxi_table.xls, xiaji_table.xls）生成run.sh, samplelist.config, new.samplelist.config
注意：目前爬表脚本对于RRBS的BIF读取可能不支持，因此RRBS需要手动填写run.sh，脚本兼容待修改。

unset PYTHONPATH
export PATH=/TJPROJ2/GB/PUBLIC/software/GB_TR/mRNA/miniconda3/bin:$PATH

python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Pipeline/pipeline/bin/bs_get_lims_info.py --name [lims登录名] --pw [密码] --pjcode [项目编号] --project_type [LiBS/RRBS/WGBS三选一]

需要修改本地已存在的final_xinxi_table.xls和xiaji_table.xls时，使用该脚本的本地版本：

python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Pipeline/pipeline/bin/bs_get_lims_info_bendi.py --name [lims登录名] --pw [密码] --pjcode [项目编号] --project_type [LiBS/RRBS/WGBS三选一]

• 微量甲基化（绝大多数情况）project_type填LiBS而非WGBS；
• RRBS注意要有–dmr_smoothing FALSE，即在DSS call DMR时不做smoothing（无生物学重复，填TRUE）；
• 生成的run.sh里第一行（ref=）需要自行填写参考基因组路径；
• –project里面的[species]和–species需要自行填写物种拉丁名；
• 需要自行填写没有爬取到的部分（BIF有可能填写不完整）；
• 需要检查确认参考基因组WGBS路径下的文件名和run.sh里的–gtf, –goann, –genenamefile, –chrlist是否对应;
• 需要检查参考基因组WGBS路径下Genome_Reg/all里是否有_utr5.bed和_utr3.bed，如果没有，则在–features里去掉utr5和utr3。

项目执行目录：/TJPROJ11/GB_TR/PJ_GB/epi/WGBS
新建一个工程文件夹，在该文件夹内按照前面的描述生成run.sh脚本，并生成samplelist.config文件，然后执行run.sh脚本：

sh run.sh

主要作用：
（1）参数整理、检验（包括文件存在性，参数有效性以及默认参数处理）；
（2）输出参数表project.config，用以检查参数；
（3）生成sjm_Analysis_afterQC.sh（QC后主分析流程）。

由于甲基化样本数据测序周期长，为压缩分析时间，目前甲基化QC新流程使用下机文库分批分lane处理

cd 1.QC_Methy

Step 1: 1.QC_Methy路径下要有new.samplelist.config，格式：datapath<tab>lib_name；执行命令：

sh QC_prepare.sh

该脚本运行后会生成sjm_split_lab.sh，运行该脚本：

sh sjm_split_lab.sh

Step 2: 如果有后续文库下机，更新1.QC_Methy/new.samplelist.config文件，重复Step 1步骤，执行命令：sh 1.QC_Methy/QC_prepare.sh ，该脚本运行后会生成sjm_split_lab.sh，运行该脚本：sh sjm_split_lab.sh

<note> 数据拆分方法
【国内】
（1）新建split.config，格式同new.samplelist.config，填写本次需拆分文库信息；
（2）执行命令：sh split_rawdata_prepare.sh，自动投递脚本；
（3）将生成的splited.samplelist.config追加至new.samplelist.config，重复step1步骤，执行命令：sh 1.QC_Methy/QC_prepare.sh ，该脚本运行后会生成sjm_split_lab.sh，运行该脚本：sh sjm_split_lab.sh。
（数据大于30G时建议拆分，且同一项目不同样本数据拆分与否需保持一致）
【海外】
样本间数据量差异过大时需要调整数据量，可根据1.QC_Methy/clean中查看样本所有份数的raw_read1文件,将不需要的某次或某几次下机的数据移除clean、raw_data备份，注意第一次下机不可移除，即*.1文件夹必须存在。移除后，sh Sample_data_summary.sh查看新数据量是否满足要求。
</note>

当所有文库都已经下机：执行命令:

sh 1.QC_Methy/QC_merged_lab.sh

该脚本运行后会生成QC_report_RUN_SJM.sh；运行该脚本：

sh QC_report_RUN_SJM.sh

QC完成后会执行QC_report.sh生成最终的QC报告。
QC完成后需检查：
<note> 1. 0.log文件夹中*_QC.JOB.status，*_QC.JOB.status.log，*_o*.txt，*_e*.txt等文件中有无报错
2. QC报告中数据量和数据质量是否达标
3. convertion_rate.txt中数字即λDNA的CT转换率高于99%
如果碰见转换率为100%的情况，首先检查bismark是否有报错，保证执行无误后检查是否为自建库，有碰见过客户未加入λDNA导致的转化率无法统计的情况。若确认客户未加入，转换率无法统计可跳过此步骤，后续分析将转换率写做100%进行分析
4. mapping率结果，若调参之后比对率仍低于60%，需向运营经理反馈
5. duplication率是否过高，内部设定阈值为30%，过高的话，若客户愿意则继续分析，因为重复reads不会进行后续分析，不会对分析造成影响，但会丢掉许多reads，导致有效数据量偏少，覆盖度难以达到要求，所以需反馈是否安排加测；客户不能接受需要重建库。
6. RRBS文库酶切效率(MspI Cutting Efficiency)高于90%
</note>

QC报告无误，上传QC报告，询问运营是否继续分析：

sh lims_QC.sh

问：BS转换率不足99%会对分析有什么影响？
答：结果假阳性会很高，不建议继续分析
问：BS转换率不足99%有什么建议？
答：重建库

在项目根目录下执行脚本：

sh sjm_Analysis_afterQC.sh

分析完成后，在项目根目录下qsub投递generate_report_result.sh，根据e、o文件检查是否有文件未能正常生成。
新流程特点：执行Analysis过程中出现有pdf但png文件为空（size=50）或不存在的情况不用管，generate_report_result.sh会统一将未生成的png由对应的pdf转换出来。

项目结题，且检查无误，qsub投递data_release.sh，生成data_release文件夹；
检查无误后，上传结题报告并释放数据：

sh lims_Report.sh

仅释放数据：

sh lims_release.sh

qsub投递backup/CMD_PJbackOSS.sh

确认结题后，qsub投递byebye.sh

格式示例：

>chr1
AGTCCTTGGATTTCCAAGA
ACACCACACCCCTTAGGAG

• 下载参考基因组及注释文件 wget *.gz
  
• “>“后为染色体名称，要与gtf文件第一列相同，且gtf第一列为fa文件chrom id的子集。
  
• 不允许存在空行。
  
• ”>“后为染色体名，无多余信息
  

chr1 ensemble exon 5843 5993 . - 0 gene_id “ENSMUSG0000171949”; transcript_id “ENSMUSG0000171949”;
chr1 ensemble CDS 5607 5806 . - 0 gene_id “ENSMUSG0000171949”; transcript_id “ENSMUSG0000171949”;

• 特点：甲基化项目gtf文件最好包括CDS, exon, intron等全部功能元件信息
  
• 作用：提取基因功能元件+DMR 结构注释+ GO+KEGG富集分析
  
• 注意：gene_id与transcript_id中经常含有”;”需要去除掉改行，否则后续会报错
  

生成方法可参考2.4（使用脚本直接生成go与genenamefile）。

ENSMUSG00000003421      GO:0019031      GO:0004970      GO:00052
ENSMUSG00000047678      GO:0004930      GO:0005515      GO:0007186      GO:0016021

• 优先从ensemble biomart下载,biomart下载不到则采用hummer2go进行注释得到
• 要求最后的GO文件不能有geneID重复的行
  

脚本及使用方法：perl /PUBLIC/source/RNA/RefRNA/version4/Prepare_scripts/delGOduplication.pl in.go out.go

• 若采用注释方法，参考/PUBLIC/source/HW/RNA/Ref/prepare_script/novelhmm.sh
  1. step1: 首先从基因组提取基因序列perl /PUBLIC/source/HW/RNA/Ref/scripts/extractcDNAfromFA.pl <gtf> <FA><OUT>
  2. step2: GO 注释：perl /PUBLIC/source/HW/RNA/Ref/scripts/hmm_pfam_go.pl -fa <基因序列文件> -n <拆分的份数>-out <输出文件前缀>
  3. step3: GO 注释结果处理,去掉行尾制表符 ：sed -i -e 's/\t$/g' *go.txt

GeneID GeneName    Description
ENSMUSG00000000001      Gnai3   Guanine nucleotide binding protein (G-protein), alpha subunit
ENSMUSG00000000003      Pbsn    Calycin-like||Odour-binding protein
ENSMUSG00000000028      Cdc45   CDC45 family

• 包含3列信息(需要含有文件头)： gene ID , gene name, gene Description

可使用以下脚本直接生成go与genenamefile。

python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Module/0.Prepare/hisat2index_go_genenamefile.py \
        --fa /TJPROJ6/GB_TR/reference_data/Plant/Osmanthus_fragrans_Lour/Osmanthus_fragrans_Lour.fa \
        --gtf /TJPROJ6/GB_TR/reference_data/Plant/Osmanthus_fragrans_Lour/Osmanthus_fragrans_Lour.gtf \
        --dir /TJPROJ6/GB_TR/reference_data/Plant/Osmanthus_fragrans_Lour/WGBS \
        --n 2,3

genenamefile生成脚本未能提取gene names,需执行以下脚本生成最终可使用的genenamefile.txt:

python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Module/0.Prepare/genenamefile.py

染色体名称，每个染色体为一行：

cat genome.gtf | cut -f 1 | uniq > chrlist.txt

生成方式例如：

cat genome.fa.fai | sort -n -k 2 | tac | head -20 | cut -f 1 > chrlist.txt

HW:/ifs/TJPROJ3/HW/PUBLIC/source/RNA/WGBS_V7/Module/0.Prepare/Ref_prepare_v2.2.py
GN:/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_RRBS_GN/WGBS/wgbs_v7.6/Module/0.Prepare/Ref_prepare_v2.3.py
统一版：/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_GB/wgbs_v7.6/Module/0.Prepare/Ref_prepare_v2.3.py
最新整合版：/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Module/0.Prepare/Ref_prepare_v2.3.py

python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Module/0.Prepare/Ref_prepare_v2.3.py \
        --genome $PWD \
        --gtf $PWD/genome.gtf \
        --goann $PWD/go.txt \
        --format gtf \
        --CGI yes \
        --repeats yes \
        --species Mus_musculus

执行完此步后，会生成sjm投递脚本,执行任务包括：

• 1. 使用bismark 软件为参考基因组建立索引；
• 2. 根据所提供的基因组注释文件，从中提取各种功能元件的信息，生成的*bed文件位于项目路径下的Genome_Reg子文件夹中；
• 3. 生成sjm_prepare.sh脚本, Genome_Reg和0.log文件夹；
• 4. 根据具体的项目要求，决定是否进行CGI 和repeats预测。

投递方法：

sh sjm_prepare.sh

将消耗内存大的任务投递到小节点，保留原来执行中任务，更换输出目录将任务投递到sm512或者smp1024节点。
/opt/gridengine/bin/linux-x64/qsub -V -cwd -l vf=*G,p=* -P smp512 *.sh

如果有个别脚本无法成功投递，投递队列中出现-q rna.q -q all.q -q novo.q -q hwus1.q -1 hwus2.q，
更改为 -q gbtr1s.q -q gball1s.q -q gball1m.q 重新投递。（需向负责人反映该情况）

如果个别脚本由于环境变量问题跑断，可以在脚本里
source /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_V7/bash_profile，重新投递

原因：RRBS文库片段较小，read中可能含有大量adapter，所以需要对adapter进行trim。有时候由于adapter.fa内缺少read1的内容，
导致read1的adapter的trim不完全。并且后续mapping test及 NT 比对都是按照read1来做，出现结果异常。
解决：在adapter.fa中添加read1内容后，重新进行trim及后续分析
正常的 adapter 文件内容如下，（）处替换成index序列，如7UDI1002的index对应为 GCTCGATT
*index对应序列文件在/TJPROJ4/XJ/script/demultiplex/DB/adapter/index.fa

>Prefixp5.7_7UDI1002/2
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>PE2/2
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>Prefixp5.7_7UDI1002/1
GATCGGAAGAGCACACGTCTGAACTCCAGTCAC（GCTCGATT）ATCTCGTATGCCGTCTTCTGCTTG
>PE/1
GATCGGAAGAGCACACGTCTGAACTCCAGTCAC（GCTCGATT）ATCTCGTATGCCGTCTTCTGCTTG

Post-Bisulfite Adaptor Tagging (PBAT)是 Bisulfite 建库过程中先转化后加接头的一种建库方式。

c图的建库方式引入了random-primer，需要对reads进行trim。使用fastp；或同流程一样，使用Trimmomatic增加HEADCROP:10参数

/TJPROJ2/GB/PUBLIC/source/Seq/software/fastp --in1 S522SA_hi.1_clean_1.fq.gz --in2 S522SA_hi.1_clean_2.fq.gz --out1 S522SA_hi_trimed_1.fq.gz --out2 S522SA_hi_trimed_2.fq.gz -q 0 -u 100 -f 10 -F 10 -n 15 -l 36 -j S522SA_hi.json -h S522SA_hi.html

本例中样本数据较大，200G/sample,为节省时间对文件拆分后进行Mapping
split -l 40000000 S522SA_hi_trimed_1.fq.gz -d -a 2 S522SA_hi_splited_

##-l 按每个文件40000000行拆分，-d 生成文件后缀按数字排序，-a 后缀的序号字符长度为3，即000,001,002,003...

perl /PUBLIC/software/RNA/Bismark/bismark_0.16.3/bismark.v16.3   -q --multicore 1  -p 4   -X 700   --dovetail     --bowtie2  --prefix map  --score_min L,0,-0.2 --unmapped --pbat  \
        --path_to_bowtie /PUBLIC/software/RNA/bowtie2-2.2.5 \
        -o  /TJPROJ6/GB_TR/PJ_GB/epi/WGBS/1100/xinyao/X204SC20090528-Z01-F003.Homo_sapiens.20201013/1.QC_Methy/clean/clean_fq/S522SA_dim.000  \
        /TJPROJ6/GB_TR/reference_data/Animal/Homo_sapiens/ensembl_homo_sapiens_grch38_p12_gca_000001405_27/WGBS  -1 S522SA_dim.splited_read1_000  -2 S522SA_dim.splited_read2_000

–pbat参数说明：
This option may be used for PBAT-Seq libraries (Post-Bisulfite Adapter Tagging; Kobayashi et al., PLoS Genetics, 2012). This is essentially the exact opposite of alignments in 'directional' mode, as it will only launch two alignment threads to the CTOT and CTOB strands instead of the normal OT and OB ones. Use this option only if you are certain that your libraries were constructed following a PBAT protocol (if you don't know what PBAT-Seq is you should not specify this option). The option ‘- -bat’ works only for FastQ files and uncompressed temporary files.

perl  /PUBLIC/software/RNA/Bismark/bismark_0.16.3/bismark.v16.3   -q --multicore 1  -p 4 --non_directional   --bowtie2  --prefix SE  --score_min L,0,-0.2  \
        --path_to_bowtie  /PUBLIC/software/RNA/bowtie2-2.2.5  \
        -o  /TJPROJ6/GB_TR/SHOUHOU/wuyou/X204SC20090528-Z01-F003.Homo_sapiens.20201013/S522SA_hi/S522SA_hi_001/check_unmappingreads  \
        /TJPROJ6/GB_TR/reference_data/Animal/Homo_sapiens/ensembl_homo_sapiens_grch38_p12_gca_000001405_27/WGBS  map.S522SA_hi_1_001_unmapped_reads_1.fq.gz  map.S522SA_hi_2_001_unmapped_reads_2.fq.gz

samtools merge map.S522SA_hi_pe.bam \
S522SA_hi_000/map.S522SA_hi_1_000_bismark_bt2_pe.bam S522SA_hi_001/map.S522SA_hi_1_001_bismark_bt2_pe.bam S522SA_hi_002/map.S522SA_hi_1_002_bismark_bt2_pe.bam S522SA_hi_003/map.S522SA_hi_1_003_bismark_bt2_pe.bam S522SA_hi_004/map.S522SA_hi_1_004_bismark_bt2_pe.bam S522SA_hi_005/map.S522SA_hi_1_005_bismark_bt2_pe.bam S522SA_hi_006/map.S522SA_hi_1_006_bismark_bt2_pe.bam S522SA_hi_007/map.S522SA_hi_1_007_bismark_bt2_pe.bam S522SA_hi_008/map.S522SA_hi_1_008_bismark_bt2_pe.bam

cd 2.Map_Methy
for i in */plot_M-bias.R;do Rscript $i;done

2)CX_accumulative,violin_plot,heatmap_clustr图片不全

CX_accumulative,violin_plot(methylation_density_violin,methylation_level_violin),heatmap_clustr图片是否完整，如果有pdf但png不全，convert转换

此脚本适用于任何有pdf无png的情况，会自动check png文件并转换

cd 3.Identification_Methy
python /TJPROJ6/GB_TR/USER/fantiantong/script/conv.py --path 3.Identification_Methy

如果pdf未生成，对于CX_accumulative

for i in */summary/CX_accumulative.r;do Rscript $i;done

对于violin和heatmap，重投mC_summary_plot.sh

参考/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Pipeline/script/wgbs_split.sh（/TJPROJ11/GB_TR/USER/xialili/script/wgbs_split.sh）把cg, chg, chh 按bin拆成更小的文件，再重新分开跑dmr，然后再合并结果.

注意修改拆出的份数和n值，份数xn 需要> 最长的chromosome 坐标，否则无法包含全部数据。

首先检查temp/split_*文件夹下dmrs.cg.xls,dmrs.chh.xls,dmrs.chg.xls是否生成完整;

然后检查所有split文件夹下是否有dmr_calling_done文件;

最后投递compare1_vs_compare2_merge.sh脚本;此脚通过dmr_calling_done确认dmr是否执行成功，确认后会删除dmr_calling_done文件后进行merge；如果merge失败需要再次投递，先在每个spli文件夹下重新touch dmr_calling_done文件

国内的环境变量/TJPROJ7/EPI/personal_dir/zhangmengran/bash_profile_zmr_20220113

在项目根目录下直接执行：python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Module/9.other/genomic_feature_level_count_plot.py

Treatment_vs_Control/DM_region下有图片生成不全，注意R版本

cd DM_region/{region_all,geneupdownstream_all}
/PUBLIC/software/public/System/R-2.15.3/bin/Rscript *.r

注意DMRvenn_promoter与DMRvenn_generegion文件夹下png是否缺失

Circos_*文件夹内circos缺失的话重新投递脚本Treatment_vs_Control_circos_level.sh

DMR_circos_*文件夹内的circos缺失，重新投递脚本Treatment_vs_Control_DMR_circos.sh

circos图生成失败，原因是import cairosvg报错，需export为Python-2.7.6的python和pythonpath；此步骤流程已修改；但是如果碰到更新前刷出来的脚本，可按如下步骤生成png文件

export PATH=/PUBLIC/software/public/System/Python-2.7.6/bin:$PATH
export PYTHONPATH=/PUBLIC/software/public/System/Python-2.7.6/lib/python2.7/site-packages:$PYTHONPATH
cd DMR_circos_CHG
mkdir t-conf  L48_treatment_vs_H_treatment_CHG.circos.confg   -svg  -outputdir  temp  -outputfile  temp/L48_treatment_vs_H_treatment_CHGraw.svg    
perl /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_RRBS_GN/WGBS/wgbs_v7.6/Module/4.DM/DMR_circos/change_svg_forcircos.V2.pl  \
temp/L48_treatment_vs_H_treatment_CHGraw.svg \
/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_RRBS_GN/WGBS/wgbs_v7.6/Module/4.DM/DMR_circos/DMR_circos.biaoci.split_contxt.txt \
temp/new_L48_treatment_vs_H_treatment_CHGraw.svg CHG_TE_gene   0 0 1 41 L48_treatment_vs_H_treatment_CHG
/PUBLIC/software/public/System/Python-2.7.6/bin/python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_RRBS_GN/WGBS/wgbs_v7.6/Module/4.DM/DMR_circos/SVG_to_other_format.py  --infile  temp/new_L48_treatment_vs_H_treatment_CHGraw.svg  --format pdf  --outfile  L48_treatment_vs_H_treatment_CHG.DMR_circos_3276800.pdf

<note> 更新历史：
甲基化新版流程V7，pipeline地址更新(如下)，与旧版v6系列相比，新版主要改动有：
1.分析主题框架改变，贯彻性的分序列环境分析（CG/CHG/CHH）。
2.有无生物学重复均使用DSS软件分析DMR。
3.添加分析模块，比较组合分析，DMR分析结果注释作图。
</note>

/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_GB/Methylation_Pipeline_V7.3.2.qsubname.py

【海外】产品线1000/1100/3000

WGBS/RRBS Pipline路径：
/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_V7/Pipeline/Methylation_Pipeline.py

项目文件夹命名：
分期号.物种拉丁名（下划线连接）.启动日期

【国内】产品线2001

WGBS Pipline路径：
/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_RRBS_GN/WGBS/wgbs_v7.6/Methylation_Pipeline_V7.3.2.qsubname.py
RRBS Pipline路径：
/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_RRBS_GN/RRBS/RRBS/Pipeline/Methylation_Pipeline_V2.py

项目文件夹命名：
分期号-B分期(一般填B1)-36（部门编号）_物种拼音_项目类型（WGBS/RRBS）_启动日期

• 表1中加粗参数均为必填参数
• 做差异分析及相应功能富集可选参数：–group，–groupname,–compare,–species,–genenamefile
  
• 表2为可选参数
• –cut，截取数据时填写
• –scoremin,调整mapping参数时填写

可提供以下pdf文件：

/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/gb_epi_wgbs/Module/9.other/Method.WGBS.pdf

在run.sh中，将–dml参数填写为yes即可。

参考路径：/TJPROJ6/GB_TR/reference_data/Animal/Magicicada_septendecula/MspI_predict
Reduced Representation Bisulfite Sequencing（RRBS）是一种准确、高效、经济的DNA甲基化研究方法，通过MspI酶识别CCGG酶切位点处理DNA，富集启动子及CpG岛区域，并进行Bisulfite处理和高通量测序
非常见物种分析前需要通过电子酶切评估，来确认MspI酶处理DNA是否能够达到建库要求。方法如下：
1. 准备参考基因组fa文件
2. 配置run.sh

ref=GCA_011763675.1_ASM1176367v1_genomic.fna #参考基因组名字
dirName='Magicicada_septendecula_MspI' #新建的结果目录的名字
window=25 #检测的windowsize，通常选25
species='Magicicada_septendecula' #物种拉丁名

perl /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_V7/RRBS_MspI_predict/GBS.En.pl \
    -c /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/WGBS/WGBS_V7/RRBS_MspI_predict/MspI_step1.cfg \
    -r $ref  \
    -o $dirName \
    -i $window \
    -s $species

3. 运行run.sh，然后进入新建文件夹投递result.sh（vf=10g,p=2）
4. 结果文件在02.stat/.stat.xls（注意是.stat.xls隐藏文件，直接打开就行）
5. 下载参考路径中的Excel文件，将结果粘贴在E2：K15，修改物种，参考基因组版本，genome size等信息。提供该表时，请勿删除人的部分，生产需要把该物种和人的结果来比较，得出是否能尝试建库的结论。