[ FRiP 스코어 ] BEDTools, Samtools 이용해서 FRiP 스코어 구하기

DiffBind 를 사용하여 분석하다가, FRiP 스코어가 너무 낮아서 다른 방법을 찾아보았다. Diffbind 에서 계산해주는 FRiP는 consesnsus peak set에 dependent하기 때문에 더 낮은 score을 얻을 수 있다고 한다. 

 

FRiP

FRiP score 이란 Fraction of Reads in Peaks이다. 말 그대로, 특정 called peaks (region) 에 들어가는 reads의 비율이라고 해석 할 수있다. 이는 usable reads in significantly enriched peaks/total number of mapped reads와 같이 계산된다.

 

실제 sample의 FRiP 계산해보기

1. 맵핑된 sequencing data (BAM/SAM) 과 called peaks files (BED/narrowPeak) 이 필요하다.

 

2. Peak regions에 intersect 하는 seqeucne reads을 알아낸다. (BEDtools 혹은 SAMtools 이용가능하다)

 

3. -> 해당 peak regions에 존재하는 Reads 갯수를 구한다.

# Step 1: Install BEDTools, SAMtools (if not already installed)

# Step 2: Intersect aligned data with peaks
bedtools intersect -a aligned_data.bam -b called_peaks.bed -u > intersected_reads.bam

# Step 3: Count the number of reads in peaks
samtools view -c intersected_reads.bam

bedtools intersect 에 대해 더 설명하자면, 두개의 파일을 비교할때, default는 정확히 겹치는 부분만 extracting 하지만 여러 다른 설정을 할 수 도 있다. https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html#default-behavior 참고 바란다.

intersect — bedtools 2.31.0 documentation

# default concept 
$ cat A.bed
chr1  10  20
chr1  30  40

$ cat B.bed
chr1  15   20

$ bedtools intersect -a A.bed -b B.bed
chr1  15   20

 

4. 총 reads 갯수를 구한다. 

# Step 4: Count the total number of reads in aligned data
samtools view -c aligned_data.bam

 

5. FRiP score = # of peaks in the peak regions / # 총 reads

# Step 5: Compute FRiP score
total_reads=32242060
reads_in_peaks=2669130
FRiP_score=$(echo "scale=2; $reads_in_peaks / $total_reads * 100" | bc) # 8.00 hmm, it should be bigger than this but 

 

 

Intersected_reads.bam 예시

# example of intersected_reads.bam
K00124:539:H7NVGBBXX:7:1113:19441:60870 99      scaffold_1      1013881 44      148M    =       1014082 351     ACATACCACCCCAAACATAAGCCCTCTTGGTATCCTTT
GAGGTCCCATGGCCTCCCAGTGTTTCTATCAAAACCAGACACCTTCATTATCCACTCATTTGTCTCCTCACGGGGCCATTTTCAAGCACAGAGTTAAGAGAGTGAGGCTG    AAFFFJJJJJJJJJFJJJJJJJJJJJJJJJJJJFJF
JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ7FJJFJJJJJJJFFJJJJJJJJJJFJJJJJJJJJJJJFJJJJJFJ    MD:Z:148        PG:Z:MarkDuplicate
s     XG:i:0  NM:i:0  XM:i:0  XN:i:0  XO:i:0  AS:i:296        YS:i:300        YT:Z:CP

# 99: flag indicating a paired-end read where both reads are mapped in the forward direction 
# scaffold_1
# 1013881: leftmost mapping position 
# 44: mapping Quality score 
# 148M: CIGAR string representing - 148 based long w/o any gaps or insertions
# 1014802: ennd
# = : same as scaffold_1
# 351