Bioinformatics

collapse taxa(example genus)

qiime taxa collapse
–i-table table.qza
–i-taxonomy taxonomy.qza
–p-level 6
–o-collapsed-table table-collapsed-L6.qza

check the number of features per taxa:

qiime feature-table summarize
–i-table table-collapsed-L6.qza
–o-visualization table-collapsed-L6-summarize.qzv
–m-sample-metadata-file sample-meta.tsv


now we do diversity core-metrics (do not miss with core-metrics-phylogenetic)

qiime diversity core-metrics \ --i-table table-collapsed-L6.qza \ --p-sampling-depth 1125 \ --m-metadata-file sample-meta.tsv \ --output-dir core-metrics_no-phylogenetics

and now we can look on PCoA:

MetaGenome 분석 파이프라인에서 OTU Table을 Biom 파일 Format으로 변환하여 메모리를 줄이는 방법이 이용되고 있습니다.

Biom 포맷은 v1.0.0과 v2.1.0이 사용되고 있는데 이 안내는 v1.0.0에 대한 설명입니다.

tsv파일을 biom으로 변환하는 방법은 json형식과 hdf5형식으로 biom으로 변환이 가능합니다.

Usage

biom convert -i feature-table.txt -o feature-table.biom --table-type="OTU table" --to-json

biom convert -i feature-table.txt -o feature-table.biom --table-type="OTU table" --to-hdf5

반대의 경우. 즉 , Biom파일을 사람이 읽기 쉽게 Tsv 파일로 변환하는 방법은 다음과 같습니다.

Usage

biom convert -i feature-table.biom -o feature-table.txt --to-tsv

추가정보)

만들어진 biom 파일을 qiime2로 Import (해당 정보는 qiim2 버전에 따라 달라질 수 있음)

구 버전)

qiime tools import --input-path feature-table.biom --type 'FeatureTable[Frequency]' --source-format BIOMV100Format --output-path table.qza

확인한 바로는 --source-format BIOMV100Format은 최신 버전에서는 넣는다면 에러가 발생합니다. 따라서 다음과 같이 command를 입력하여야 합니다.

2018년 이후 버전)

qiime tools import --input-path feature-table.biom --type 'FeatureTable[Frequency]' --output-path table.qza

Chip seq 분석을 마친 후 narrowPeak 파일을 이용하여 Peak Profiling 을 진행하고 figure를 제작할 수 있습니다.

사용한 library는 ChIPseeker와 ggplot2 입니다.

narrowPeak 파일 format은 다음과 같습니다.

1. chrom - Name of the chromosome (or contig, scaffold, etc.).
2. chromStart - The starting position of the feature in the chromosome or scaffold. The first base in a chromosome is numbered 0.
3. chromEnd - The ending position of the feature in the chromosome or scaffold. The chromEnd base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined aschromStart=0, chromEnd=100, and span the bases numbered 0-99.
4. name - Name given to a region (preferably unique). Use '.' if no name is assigned.
5. score - Indicates how dark the peak will be displayed in the browser (0-1000). If all scores were '0' when the data were submitted to the DCC, the DCC assigned scores 1-1000 based on signal value. Ideally the average signalValue per base spread is between 100-1000.
6. strand - +/- to denote strand or orientation (whenever applicable). Use '.' if no orientation is assigned.
7. signalValue - Measurement of overall (usually, average) enrichment for the region.
8. pValue - Measurement of statistical significance (-log10). Use -1 if no pValue is assigned.
9. qValue - Measurement of statistical significance using false discovery rate (-log10). Use -1 if no qValue is assigned.
10. peak - Point-source called for this peak; 0-based offset from chromStart. Use -1 if no point-source called.

track type=narrowPeak visibility=3 db=hg19 name="nPk" description="ENCODE narrowPeak Example"
browser position chr1:9356000-9365000
chr1    9356548 9356648 .       0       .       182     5.0945  -1  50
chr1    9358722 9358822 .       0       .       91      4.6052  -1  40
chr1    9361082 9361182 .       0       .       182     9.2103  -1  75

아래 코드들은 R 에서 작성되었습니다.

Usage1

library(ChIPseeker)
library(ggplot2)

peak <- readPeakFile("narrowPeak파일")
peak

## GRanges object with 1331 ranges and 2 metadata columns:
##          seqnames              ranges strand |             V4        V5
##             <Rle>           <IRanges>  <Rle> |       <factor> <numeric>
##      [1]     chr1       815093-817883      * |    MACS_peak_1    295.76
##      [2]     chr1     1243288-1244338      * |    MACS_peak_2     63.19
##      [3]     chr1     2979977-2981228      * |    MACS_peak_3    100.16
##      [4]     chr1     3566182-3567876      * |    MACS_peak_4    558.89
##      [5]     chr1     3816546-3818111      * |    MACS_peak_5     57.57
##      ...      ...                 ...    ... .            ...       ...
##   [1327]     chrX 135244783-135245821      * | MACS_peak_1327     55.54
##   [1328]     chrX 139171964-139173506      * | MACS_peak_1328    270.19
##   [1329]     chrX 139583954-139586126      * | MACS_peak_1329    918.73
##   [1330]     chrX 139592002-139593238      * | MACS_peak_1330    210.88
##   [1331]     chrY   13845134-13845777      * | MACS_peak_1331     58.39
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

출처 : https://genome.ucsc.edu/FAQ/FAQformat.html#format12 

Usage2

peak <- readPeakFile("narrowPeak파일")
covplot("narrowPeak파일")

♣각각의 chr에서의 peak을 profile 할 수 있습니다.

Usage3

peak = GenomicRanges::GRangesList(샘플명=readPeakFile("narrowPeak파일"), 샘플명=readPeakFile("narrowPeak파일")
covplot(peak, weightCol="V5") 

♣Control과 Test의 peak의 차이를 한번에 확인할 수 있습니다.

Usage4

peak = GenomicRanges::GRangesList(샘플명=readPeakFile("narrowPeak파일"), 샘플명=readPeakFile("narrowPeak파일")
covplot(peak, weightCol="V5", chrs = c("chr선택","chr선택"), xlim=c(x축 start position, X축 end position)) 

♣ 원하는 각 chromosome의 size를 제한 하여 peak을 볼 수 있습니다.

Additional Options

xlab  -> x축 title

ylab -> y축 title

title -> figure title

----------------------------------

### + facet_grid(chr ~.id)  -> 범례를 설정합니다. 또한 y축의 scale이 setting 됩니다.