Revision as of 04:27, 22 June 2020

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

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.vcf
samtools 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은 아마도 업그레이드 된 형태일 듯

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

분석 과정:

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 만들기, 이후 모부본 자리 바꾸기

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