1: Reference 구축
*해당 문서는 CellRanger-ATAC/1.2.0 기준으로 작성되었습니다. (최신 업데이트 버전 :2.0)
1. fasta file과 gtf 파일 다운로드
- Ensembl: 모든 종
- phytozomev11: 식물의 경우(Maize는 maizegdb.org)
2. Config 파일 만들기
CellRanger 는 처음에 fasta file 과 gtf파일을 다 받아서 매핑을 하고 뒤에 매핑된 리드들이 어디 부분에 위치되어있는지까지 알려주는 파이프라인이기 때문에,
처음 매핑할 때 여러 정보가 요구된다.
2.0 버전의 경우 config파일의 형식이 약간 다를 수 있으며, gtf파일에서 exon정보 및 transcript_id등을 요구한다. (https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/advanced/references#jaspar)
따라서 해당 gtf가 없는 종을 연구하고 있다면 1.2.0을 사용하는게 무난하다.
vi명령어를 이용하여 config 파일하나를 만들어 준다.
<1.2 ver>
{ GENOME_FASTA_INPUT: "파일경로/Athaliana_447_TAIR10.fa", GENE_ANNOTATION_INPUT: "파일경로/Athaliana_447_Araport11.gene_exons.gff3", MOTIF_INPUT: "", ORGANISM: "Arabidopsis", PRIMARY_CONTIGS: ["Chr1", "Chr2", "Chr3", "Chr4", "Chr5"],
NON_NUCLEAR_CONTIGS: ["ChrM","ChrC"]}
<2 ver>
{ organism: "Populus"
genome: ["populus"]
input_fasta: ["/0.Reference_populus/g717v5_h1h2.fa"]
input_gtf: ["/0.Reference_populus/PtremulaxPopulusalbaHAP1v5.1.gene_exons.gtf"]}
만약에 가지고있는 어노테이션 파일이 gff3 라면 마지막 attribute 열의 이름이 일치하지 않아 에러가 뜰수 있다.
이때는 이전 글 "gff3 --> gtf 변환 포스트"에서 gtf로 변환방법을 통해 gtf로 변환을 먼저 하길 추천한다.
gtf 파일 상태에 따라서 ver2를 써야하는 경우도 있는데 이때 config 파일이 다르다.
3. Reference 구축
<1.2 ver>
module load CellRanger-ATAC/1.2.0 cellranger-atac mkref 아웃풋파일경로(꼭새로운폴더여야함!) --config=내가 만든 config파일 이름
소요시간: 1시간
<2 ver>
module load CellRanger-ATAC/2.0.0 cellranger-atac mkref --config=내가 만든 config파일 이름
만약에 샘플이 Haploid면 어떻게 할까?
2: Read mapping
1. Single Cell RNA-seq/ATAC-seq data 이해
일반 시퀀싱 파일은 single end로 시퀀싱을 할경우 하나의 ".fastq"파일을 갖게되고,
paired end 의 경우 "R1.fastq","R2.fastq"두개의 파일을 갖게됩니다.
Single cell의 경우 4개의 파일을 갖게되고 이름은 각각
"I1.fastq","R1.fastq","R2.fastq" "R3.fastq"입니다.
"I1": 8bp길이의 바코드
"R1": Forward end read
"R2": 16bp길이의 10x feature바코드
"R3": Reverse end read
참고: https://divingintogeneticsandgenomics.rbind.io/post/understand-10x-scrnaseq-and-scatac-fastqs/
2. cellranger-atac count를 통한 mapping
소요시간: 32코어를쓸 때, 20gb fastq 20시간
<Ver 1.2. 와 Ver 2 동일>
cellranger-atac count \ --id=ArabRoot-rep1-scATAC \ --reference=1단계에서 구축한 reference경로 \ --fastqs=파스타파일이 있는경로 --localcores=2
3: 중복 맵핑된 reads를 제거
중복 맵핑된 Reads를 제거
어떤 Read들은 2개이상의 genome 부분에 붙을 수 있다. 이러한 Read들을 제거해 주는 단계이다.
samtools view -@ 32 -h -f 3 -q 20 ./"${List[SLURM_ARRAY_TASK_ID]}"/outs/possorted_bam.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -@ 32 -bS - > ../3.SortedBam/"${List[SLURM_ARRAY_TASK_ID]}"_Sorted.bam
옵션 설명
-h: output에 헤더 포함
-f 3: flag중 3만 출력. https://broadinstitute.github.io/picard/explain-flags.html.
3은 paired read 중 적절하게 맵핑된 Read들만을 이야기한다.
-q 20: Skip alignments with MAPQ smaller than INT
4: PCR duplicates 제거
java -jar $EBROOTPICARD/picard.jar MarkDuplicates \ MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 \ REMOVE_DUPLICATES=true \ METRICS_FILE= ./"${List[SLURM_ARRAY_TASK_ID]}"_dups.txt(아웃풋) \ I= ./"${List[SLURM_ARRAY_TASK_ID]}"_Sorted.bam(이전단계-중복제거된 뱀파일) \ O= ./"${List[SLURM_ARRAY_TASK_ID]}"_Rmpcr.bam(아웃풋) \ BARCODE_TAG=CB \ ASSUME_SORT_ORDER=coordinate\ USE_JDK_DEFLATER=true USE_JDK_INFLATER=true
이때 만약에 앞에서 mapping 을 위해 cell ranger버젼 2를 사용한 경우엔 header 에러가 뜨게 된다.
이때는 samtools view -H 를 이용해서 bam file 의 헤더를 관찰할 수 있다.
samtools view -H /3.scATAC_flo/3.SortedBam/Populus_Sorted.bam > Populus_Sorted_header.sam
@hd SO:coordinate @sq SN:T01 LN:50695696 @sq SN:T02 LN:23110994 @sq SN:T03 LN:25238799 @sq SN:T04 LN:21613846 @sq SN:T05 LN:22266242 @sq SN:T06 LN:25485240
이때, 첫번째 헤더의 경우 VN version number 가 빠져있다. 이때는 먼저
위의 헤더의 "Populus_Sorted_header.sam"를 열고 안에 내용을 바꾸어준다.
vi Populus_Sorted_header.sam
@hd VN:1.5 SO:coordinate
@sq SN:T01 LN:50695696 @sq SN:T02 LN:23110994 @sq SN:T03 LN:25238799 @sq SN:T04 LN:21613846
그 다음에 아래 명령어를 이용해서 bamfile의 헤더를 바꾼다.
samtools reheader /3.scATAC_flo/3.SortedBam/Populus_Sorted_header.sam /3.scATAC_flo/3.SortedBam/Populus_Sorted.bam > /3.scATAC_flo/3.SortedBam/Populus_Sorted_HF.bam
5: Barcode 고정
6: Unique Read만 bed file 로 가져오기
1) makeTn5bed.pl 를이용!
단, 이때 해당 펄 파일은 perlSAM.pm 라는 모듈을 필요로 하기 때문에,
모듈과 펄 실행파일을 같은 디렉토리에 넣고,
conda를 통해 펄을 실행시킨다 (확인은 which perl)
그리고 Naturally라는 펄 모듈은 콘다로 다운받아줘야한다.
IDR scores
LCM data sets
