CB population genetic map
WD: 244:/hayasen/chojam/new_GBS
하위 디렉토리:
1) GBS --> /alima9002/YMY/UV/GBS
태영이형이 윤박사님 데이터 분석한 폴더
176 line들 중 138라인만 .bam.sorted.bam 형태로 되어 있고,
나머지 38라인과 모,부본은 .cat.fq.trimmed 형태로 되어 있다.
2) UV --> /NGS/NGS/GlycineMax/DNA/UV
윤박사님이 DNA추출 후 얻은 GBS raw data
분석은 barcode removal까지 된 .fq_trimmed 형식으로 된 파일들을 사용한다.
3) assembly --> /alima9002/ref/Gmax/assembly
Gmax275 genome sequence
Contents |
1st Trial
분석 과정:
1) the_38_lines.list는 모종의 이유로 정상적인 분석에서 제외된 38라인과 모,부본을 담고 있다
이번 분석에서는 일단 이 라인들을 제외하지 않고, 다시 reference에 mapping을 하여 bamfile을 생성했다.
예시:
bwa mem (-t 4) Gmax_275_v2.0.fa 6_barcode.lane7_UV77.fq_trimmed > UV77.bam samtools sort -o UV77.bam.sorted.bam UV77.bam
2) 그렇게 총 174 + 모부본의 bam.sorted.bam 파일이 생성되고, mpileup을 진행한다.
(원래는 samtools mpileup을 썼던 것 같은데 이제는 bcftools mpileup을 사용하나보다)
command:
bcftools mpileup -b bam_list -Ou -f Gmax_275_v2.0.fa | bcftools call -Ov -mv -o UV.chojam.variant.vcfsamtools mpileup -f Gmax_275_v2.0.fa -v -t DP,AD,ADF,ADR,SP,INFO/AD,INFO/ADF,INFO/ADR -u -b bam_list | bcftools call -Ov -mv -o UV.chojam.variant.vcf
samtools mpileup의 -t 옵션은 bcftools mpileup의 -a (annotate) 옵션으로 치환가능하다.
*options: -Ou: 중간단계에서 textual 형태로 쓰지 않고, uncompressed binary 형태로 가지고 있어라 -Ov: vcf 형태의 파일을 만들어라 -v: variant만 calling해라 -m: -c 옵션이 예전 mpileup 옵션이었다면, -m은 아마도 업그레이드 된 형태일 듯
-> -O, --output-type b|u|z|v Output compressed BCF (b), uncompressed BCF (u), compressed VCF (z), uncompressed VCF (v). Use the -Ou option when piping between bcftools subcommands to speed up performance by removing unnecessary compression/decompression and VCF←→BCF conversion.
3) 생성된 vcf 파일 filtering - 일단은 depth, quality에 대해서만 진행
command:
vcftools --vcf UV.chojam.variant.vcf --out UV.chojam.d3.q30 --remove-indels --minDP 3 --minQ 30 --recode
4) 모부본에 대해서 homo이고 polymorphic한 SNP만 걸러낸다.
command:
python parent_homo.py UV.chojam.d3.q30.recode.vcf UV.chojam.d3.q30.recode.homo.vcf
5) missing option: max-missing 0.1 0.3 0.5 0.7
command:
vcftools --vcf UV.chojam.d3.q30.recode.homo.vcf --out UV.chojam.d3.q30.m0.1 --recode --max-missing 0.1
--max-missing 0.1 -> 10178 --max-missing 0.3 -> 9277 --max-missing 0.5 -> 7599 --max-missing 0.7 -> 4899
6) transform into .loc file
command:
python vcf_to_loc.py [vcffile] [population_name] [population_type] [out_locfile] python vcf_to_loc.py UV.chojam.d3.q30.m0.5.recode.vcf new_CB RI6 UV.chojam.d3.q30.m0.5.recode.vcf.loc
근데 이걸 보고 있으니, 몇 개의 라인에서 계속적으로 genotype이 안 읽힌다. 아마 '그' 38개일 듯
2nd Trial
Data processing
분석 과정:
1) '그' 38개 파일의 출처 확인
리드 수가 안 나와서 두 번에 걸쳐 GBS 한 것으로 판단됨
첫 번째 GBS는 '/NGS/NGS/GlycineMax/DNA/UV'에 초록색으로 표시된 '#_barcode.lane#_UV##.fq_trimmed' 형태의 파일이고,
두 번째 GBS ('the_38_lines')는 같은 디렉토리에 있는 하얀색 '#_UV##A.fq_trimmed' 형태의 파일이다.
'그' 38개의 라인은 첫 번째와 두 번째 파일이 합쳐진 'UV##.cat.fq.trimmed' 형태이다.
2) 태영이형이 이미 mapping까지 끝낸 파일들이 있으므로, 그것들을 사용한다.
command:
cp /alima9002/YMY/UV/GBS/*.bam.sorted.bam ./new_bamfiles -> 모부본 포함 총 176개 파일 ls -1 | sort -V > new_bam_list -> 새로운 bam_list 만들기, 이후 모부본 자리 바꾸기
*option: sort -V: sort by lexicographic order(!)
3) mpileup
command:
bcftools mpileup -b new_bam_list -a DP,AD,ADF,ADR,SP,INFO/AD,INFO/ADF,INFO/ADR -Ou -f Gmax_275_v2.0.fa | bcftools call -Ov -mv -o UV.new.chojam.variant.vcf
4) filtering -> 두 개의 파일 생성: InDel 없는 파일 (SNP-only), InDel만 있는 파일
*InDel의 경우 quality를 조금 높게 (>=50) 설정하여 integrity를 높인다 *depth도 높게 (>=5) 설정한다
command:
vcftools --vcf UV.new.chojam.variant.vcf --out UV.new.chojam.d3.q30.SNP --minDP 3 --minQ 30 --remove-indels --recode vcftools --vcf UV.new.chojam.variant.vcf --out UV.new.chojam.d5.q50.indel --minDP 5 --minQ 50 --keep-only-indels --recode
5) run parent_homo.py
command:
python parent_homo.py UV.new.chojam.d3.q30.SNP.recode.vcf UV.new.chojam.d3.q30.SNP.recode.homo.vcf python parent_homo.py UV.new.chojam.d5.q50.indel.recode.vcf UV.new.chojam.d5.q50.indel.recode.homo.vcf
6) missing option: max-missing 0.2 low stringency 0.5 high stringency
command:
vcftools --vcf UV.new.chojam.d3.q30.SNP.recode.homo.vcf --out UV.new.chojam.d3.q30.m0.2.SNP --recode --max-missing 0.2 (0.5) vcftools --vcf UV.new.chojam.d5.q50.indel.recode.homo.vcf --out UV.new.chojam.d5.q50.m0.2.indel --recode --max-missing 0.2 (0.5)
7) 최종 파일은 SNP과 InDel을 합친 파일
m = 0.2 -> UV.new.chojam.d3.q30.m0.2.SNP.recode.vcf + UV.new.chojam.d5.q50.m0.2.indel.recode.vcf -> UV.new.chojam.m0.2.SNPwINDEL.vcf (12469)
m = 0.5 -> UV.new.chojam.d3.q30.m0.5.SNP.recode.vcf + UV.new.chojam.d5.q50.m0.5.indel.recode.vcf -> UV.new.chojam.m0.5.SNPwINDEL.vcf (9657)
8) transform into .loc file
command:
python vcf_to_loc.py [vcffile] [population_name] [population_type] [out_locfile] python vcf_to_loc.py UV.new.chojam.m0.2.SNPwINDEL.vcf new_CB RI6 UV.new.chojam.m0.2.SNPwINDEL.vcf.loc python vcf_to_loc.py UV.new.chojam.m0.5.SNPwINDEL.vcf new_CB RI6 UV.new.chojam.m0.5.SNPwINDEL.vcf.loc
Linkage map construction
0) 오류 감지
UV.new.chojam.m0.5.SNPwINDEL.vcf.loc 이 파일을 JoinMap에 옮겼더니 중복된 마커가 감지됨.
double presence of locus Chr08_46691793 (numbers 3193 and 9244)
1) 그래서 일단은 각각의 SNP, InDel 마커를 구분한다. vcf_to_loc.py를 다듬는다.
command:
python vcf_to_loc.py UV.new.chojam.m0.2.SNPwINDEL.vcf new_CB RI6 UV.new.chojam.m0.2.SNPwINDEL.vcf.loc python vcf_to_loc.py UV.new.chojam.m0.5.SNPwINDEL.vcf new_CB RI6 UV.new.chojam.m0.5.SNPwINDEL.vcf.loc
2) JoinMap 옵션들
일단 이하 JoinMap 과정에 들어가는 파일은 missing 50%인 마커들이다.
- 2-1) Similarity of Loci -> Population -> Exclude identicals 1041 identical loci excluded + Exclude similar loci (0.995) No loci excluded
- 2-2) Locus Genotype Frequency -> Highly significant loci excluded (prone to segregation distortion) 유의성 4개 (****) 이상인 것 제외
JoinMap grouping 중, LOD 25.0로 설정했으나 갈라지지 않은 마커가 발견됨 -> vcf filtering을 더 빡세게 해야겠다고 생각함
vcf filtering options: QUAL >= 999, DP >= 5, max-missing = 0.2/0.5
vcftools --vcf UV.new.chojam.m0.2.SNPwINDEL.vcf --out UV.new.chojam.d5.q999.m0.2.SNPwINDEL --minDP 5 --minQ 999 --max-missing 0.2 --recode vcftools --vcf UV.new.chojam.m0.5.SNPwINDEL.vcf --out UV.new.chojam.d5.q999.m0.5.SNPwINDEL --minDP 5 --minQ 999 --max-missing 0.5 --recode
JoinMap options: similarity >= 0.95, segregation distortion >= ***
3) 총 24개의 group으로 grouping, prefix snpwindel_m50_d5q999로 저장 (JoinMap, Excel) -> ML interval mapping
