m6A-MeRIP and m6Am-Exo-MeRIP-seq analysis

Fastqc was used to perform quality control on sequencing data. And then, cutadapt was applied to remove adapters and trim reads. The pre-processed reads were aligned to human genome (built on Release 31 (GRCh38.p12) gencode) by STAR⁵³. To call m⁶A peaks, MACS2⁵⁴ (q value: 0.05, call-summit mode) was used based on its paired m⁶A-RIP/input data. The default peak range of m⁶A viewer is 200 nt and above 80% of peaks form MACS2 has range around 200 nt. To find the collection of common peaks for each group, peaks from biological replicates are filtered using an adjusted method based on the previously published paper⁵⁵. Briefly, to define the common peaks, they should appear in at least two biological replicates, and peaks within 200 nt from each other are defined as one peak. These kept peaks together build the collection of common peaks for each group. To generate the peak distribution across chromosome region (including intron, CDS, intergenic region) and across gene region (3′UTR, 5′UTR and CDS) RSeQC⁵⁶ and Guitar Plot⁵⁷ were applied. The gene annotation used here is from GENCODE⁵⁸. For further m⁶Am peak analysis, we extracted common peaks whose summit is in the 5′UTR region. The data visualization of these peaks is realized using bedtools⁵⁹ and IGV^60,61. To find potential m⁶Am peaks, peaks in 5′UTR region in the control group were compared with PCIF1 KO groups. In MeRIP seq, peaks in controlled groups but not in two PCIF1 KO groups (HIV− and HIV+) were considered as potential m⁶Am peaks. These potential m⁶Am peaks were then compared with peaks in HIV-infected groups to identify HIV-changed peaks. Potential peaks which have decreased peaks in control and HIV-infected groups are considered as m⁶Am peaks changed by HIV infection.