21： Jing-ru EM-seq 解析 (Lib1〜8)

###### tags:`ノート` `JL` `解析` `EM-seq` 21： Jing-ru EM-seq 解析 (Lib1〜8) === #### 【作業場所＆保存場所】 //osc-fs_home/s-maeda/21_JL_EM-seq_20220316 #### 【rawdata】 //analysisdata/rawseq/fastq/SHARED/000037/21_JL_EM-seq_20220316_rawdata #### 【参考】 bismarkでのメチル化解析 HackMD https://hackmd.io/@takahirosuzuki/B1VFrdwoH --- #### 2022/03/16 ### 1. Adapter Trimming 今回の配列 Read 1: AGATCGGAAGAGCACACGTCTGAACTCCAGTCA Read 2: AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT #### 1_cutadapt.sh ``` #!/bin/bash output1=$1 output2=$2 read1=$3 read2=$4 cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT -o $output1 -p $output2 $read1 $read2 ``` ### 2. クオリティコントロール #### 2_prinseqpp.sh ``` #!/bin/bash read1=$1 read2=$2 output=$3 conda activate prinseq-plus-plus prinseq++ -fastq $read1 -fastq2 $read2 -out_name $output -trim_tail_left 5 -trim_tail_right 5 -trim_qual_left 30 -trim_qual_right 30 -ns_max_n 20 -min_len 50 -threads 8 ``` ### 3. FastQC #### 3_fastqc.sh ``` #!/bin/bash read=$1 fastqc --nogroup -t 8 -o FastQC_reports $read ``` FastQCの結果は //osc-fs_home/s-maeda/21_JL_EM-seq_20220316/FastQC_reports にある。 ### 4. Bismark マッピング #### 4_bismark.sh ``` #!/bin/bash read1=$1 read2=$2 bismark --bowtie2 --fastq --parallel 8 -N 1 bisulfite -1 $read1 -2 $read2 ``` Lib1 Mapping efficiency: 70.8% Lib2 Mapping efficiency: 53.0% Lib3 Mapping efficiency: 72.4% Lib4 Mapping efficiency: 74.0% Lib5 Mapping efficiency: 76.9% Lib6 Mapping efficiency: 76.9% Lib7 Mapping efficiency: 76.8% Lib8 Mapping efficiency: 76.4% Lib2のマッピング効率が低い。アダプター少し変えて再度やってみたが変わらず... ### 5. duplicateの除去 ``` #!/bin/bash bamfile=$1 deduplicate_bismark --bam -p $bamfile ``` ### 6. Methylation call ``` #!/bin/bash bamfile=$1 bismark_methylation_extractor --gzip --bedGraph --parallel 4 --buffer_size 10G --cytosine_report --genome_folder //osc-fs_home/s-maeda/21_JL_EM-seq_20220316/bisulfite $bamfile ``` ### 7. IGVで可視化用の材料つくるソートしてindexを作る ``` #!/bin/bash #SBATCH -p amd #SBATCH -n 4 bamfile=$1 outputbam=$2 samtools sort $bamfile -o $outputbam samtools index $outputbam ``` *.bedGraph, ソートしindexを作ったbamファイルをIGVなどのgenome browserで読み込むことで可視化できる。 ### 8. MethylKitでそれぞれ比較 #### Sample情報 | Number | Name | | ---- | ---- | | Lib1 | WT_Day6　　| | Lib2 | WT_Day9 | | Lib3 | KO_Day6 | | Lib4 | KO_Day9 | | Lib5 | VP64_Day6 | | Lib6 | VP64_Day6 | | Lib7 | WT_Day9_P6 | | Lib8 | VP64_Day9_P6 | #### 組み合わせ * WT_Day9　vs　WT_Day6 * KO_Day9　vs　KO_Day6 * VP64_Day9　vs　VP64_Day6 * WT_Day9_P6　vs　WT_Day6 * VP64_Day9_P6　vs　VP64_Day6 ### bam fileの読み込み 8_methylKit_read2.sh ``` #!/bin/bash #SBATCH -p amd #SBATCH -n 4 source ~/.bash_profile R --vanilla --slave --args $1 $2 $3 $4 $5 < 8_methylKit_read2.R ``` 8_methylKit_read2.R ``` args <- commandArgs(trailingOnly = T) library(methylKit) file1 <- args[1] file2 <- args[2] sample1 <- args[3] sample2 <- args[4] rdata <- args[5] files <- list(file1, file2) sample.id <- list(sample1, sample2) group_vct <- c(0,1) my.methRaw.list=processBismarkAln( location = files, sample.id=sample.id, assembly="hg38", read.context="CpG", save.folder=getwd(), mincov = 4, treatment=group_vct) save(list=ls(), file=rdata) ```