Difference between revisions of "Methylation study"
KangHeum Cho (Talk | contribs) (→HTSeq) |
KangHeum Cho (Talk | contribs) |
||
| Line 1: | Line 1: | ||
= Tuxedo = | = Tuxedo = | ||
| + | |||
| + | WD: 63:/data6/chojam96/methylation/tuxedo | ||
| + | |||
== Basic Bowtie2 pipeline == | == Basic Bowtie2 pipeline == | ||
| Line 26: | Line 29: | ||
* ''Name'' for GFF format | * ''Name'' for GFF format | ||
* etc | * etc | ||
| + | |||
| + | == edgeR == | ||
| + | |||
| + | Refer to [[2108 Hakyung TUXEDO]] | ||
| + | |||
| + | If edgeR not installed properly due to some version issues, follow instructions from [https://bioconductor.org/packages/release/bioc/html/edgeR.html edgeR installation]. | ||
| + | |||
| + | |||
| + | = DEG and methylation = | ||
| + | |||
| + | == DEG == | ||
| + | |||
| + | From HTSeq, for six samples (three for K, three for S), files ''*_htseq.txt'' were produced and combined into ''combined_htseq.txt'' | ||
| + | |||
| + | edgeR was performed on 'combined_htseq.txt' on R studio and 'combined_htseq.count.edgeR' (n=23106) was produced. | ||
| + | |||
| + | DEG was extracted using: | ||
| + | python finddeg.py combined_htseq.count.edgeR combined_htseq.deg | ||
| + | |||
| + | DEGs from each cytosine contexts were extracted using: | ||
| + | python extractmetccount.py <DEG file> <*_genes> <*_count.deg> | ||
| + | |||
| + | CpG_count.deg (n=21) | ||
| + | |||
| + | CHG_count.deg (n=23) | ||
| + | |||
| + | CHH_count.deg (n=4) | ||
| + | |||
| + | == Gene methylation count == | ||
| + | |||
| + | WD: 63:/data6/chojam96/methylation/expression_methylation | ||
| + | |||
| + | |||
| + | 1. From ''193:/data2/chojam96/methylation/metilene/metilene_v0.2-7/Vr-Met-S(K)'', get bismark coverage files | ||
| + | |||
| + | Six files: 'CpG(H)_Vr-Met-S(K)_bismark.cov' | ||
| + | |||
| + | Chr Pos Pos MetPercentage Met_read Unmet_read | ||
| + | Vr01 44 44 100 1 0 | ||
| + | Vr01 45 45 100 1 0 | ||
| + | Vr01 359 359 100 1 0 | ||
| + | Vr01 360 360 100 2 0 | ||
| + | Vr01 371 371 100 1 0 | ||
| + | |||
| + | 2. Combine all cytosine context files → 'Vr-Met-S(K)_bismark.cov' and set minimum threshold as (met + unmet) >= 3 → 'Vr-Met-S(K)_bismark.min3.cov' | ||
| + | |||
| + | Chr Pos Met_read Unmet_read | ||
| + | Vr01 417 4 0 | ||
| + | Vr01 459 2 1 | ||
| + | Vr01 467 1 2 | ||
| + | Vr01 606 6 0 | ||
| + | Vr01 608 3 1 | ||
| + | |||
| + | 3. Divide into chromosomes and scaffolds → 'Vr#(scaffold)_S(K).min3.cov' | ||
| + | |||
| + | python dividescaffolds.py Vr-Met-S(K)_bismark.min3.cov | ||
| + | |||
| + | 4. '''Run 'getdata.py'''': takes GFF - only mRNA - file, complete gene ID and line type (S, K) | ||
| + | |||
| + | '''Coverage(.cov) files of ''all the scaffolds'' should be present in the same working directory!''' | ||
| + | |||
| + | python getdata.py <GFF.onlymRNA> <Gene ID> <nowiki><S/K></nowiki> | ||
| + | python getdata.py ../Vradi_ver6.onlymrna.gff Vradi01g05390.1 S | ||
| + | |||
| + | This produces <Gene ID>_<nowiki><S/K></nowiki> file with following format: | ||
| + | |||
| + | 'Vradi01g05390.1_S' | ||
| + | ''Position Number of counts Met/Unmet'' | ||
| + | #strand:- | ||
| + | 8291109 4 met | ||
| + | 8291109 -0 unmet | ||
| + | 8291141 5 met | ||
| + | 8291141 -0 unmet | ||
| + | 8291162 4 met | ||
| + | 8291162 -2 unmet | ||
| + | 8291182 12 met | ||
| + | 8291182 -0 unmet | ||
Revision as of 00:46, 25 April 2019
Contents |
Tuxedo
WD: 63:/data6/chojam96/methylation/tuxedo
Basic Bowtie2 pipeline
Refer to 2108 Hakyung TUXEDO
HTSeq
For detailed information, refer to HTSeq page.
Basic command:
htseq-count -t exon (-f bam) (--stranded=no) (-q) <sorted.bam> <annotation.gff/gtf> > <htseq.txt>
Default input file format is SAM, so specifying '-f bam' for BAM input file is necessary.
For files not following normal HTSeq pipeline, '--stranded=no' option is also required. Refer to 2108 Hakyung TUXEDO.
However, for V. radiata Illumina-sequenced GFF3 files, ID attributes are mixed up, so additional '--idattr' option should be used!
htseq-count -t exon -f bam --stranded=no -q --idattr transcript_id <sorted.bam> <annotation.gff/gtf> > <htseq.txt>
ID attributes may be:
- transcript_id for GTF format
- Name for GFF format
- etc
edgeR
Refer to 2108 Hakyung TUXEDO
If edgeR not installed properly due to some version issues, follow instructions from edgeR installation.
DEG and methylation
DEG
From HTSeq, for six samples (three for K, three for S), files *_htseq.txt were produced and combined into combined_htseq.txt
edgeR was performed on 'combined_htseq.txt' on R studio and 'combined_htseq.count.edgeR' (n=23106) was produced.
DEG was extracted using:
python finddeg.py combined_htseq.count.edgeR combined_htseq.deg
DEGs from each cytosine contexts were extracted using:
python extractmetccount.py <DEG file> <*_genes> <*_count.deg>
CpG_count.deg (n=21)
CHG_count.deg (n=23)
CHH_count.deg (n=4)
Gene methylation count
WD: 63:/data6/chojam96/methylation/expression_methylation
1. From 193:/data2/chojam96/methylation/metilene/metilene_v0.2-7/Vr-Met-S(K), get bismark coverage files
Six files: 'CpG(H)_Vr-Met-S(K)_bismark.cov'
Chr Pos Pos MetPercentage Met_read Unmet_read Vr01 44 44 100 1 0 Vr01 45 45 100 1 0 Vr01 359 359 100 1 0 Vr01 360 360 100 2 0 Vr01 371 371 100 1 0
2. Combine all cytosine context files → 'Vr-Met-S(K)_bismark.cov' and set minimum threshold as (met + unmet) >= 3 → 'Vr-Met-S(K)_bismark.min3.cov'
Chr Pos Met_read Unmet_read Vr01 417 4 0 Vr01 459 2 1 Vr01 467 1 2 Vr01 606 6 0 Vr01 608 3 1
3. Divide into chromosomes and scaffolds → 'Vr#(scaffold)_S(K).min3.cov'
python dividescaffolds.py Vr-Met-S(K)_bismark.min3.cov
4. Run 'getdata.py': takes GFF - only mRNA - file, complete gene ID and line type (S, K)
Coverage(.cov) files of all the scaffolds should be present in the same working directory!
python getdata.py <GFF.onlymRNA> <Gene ID> <S/K> python getdata.py ../Vradi_ver6.onlymrna.gff Vradi01g05390.1 S
This produces <Gene ID>_<S/K> file with following format:
'Vradi01g05390.1_S' Position Number of counts Met/Unmet #strand:- 8291109 4 met 8291109 -0 unmet 8291141 5 met 8291141 -0 unmet 8291162 4 met 8291162 -2 unmet 8291182 12 met 8291182 -0 unmet
