BSA分析:GATK 4の使用(bwaを含む)
21045 ワード
注:バージョンが異なると、コマンドが一致しません.必ず対応するバージョンを使います.
GoogleアカウントGATK 3のダウンロードアドレスGATK 4の各バージョンにログインする必要があり、コマンドによっては変更があります.1.GATKのインストール、他人のチュートリアルテンセントクラウドを使用するチュートリアルGATK 4データ前処理変異検出(BWA+SAMtools+picard+GATK)bwa+samtools+picardtools+GATK call SNPプロセスは、私たちのサーバーのスラグのネットワーク問題のため、内容はすべて地元のwin 10にダウンロードした後、サーバーにアップロードします.サーバ構成:cpu 4個6コア2スレッド=48メモリ:396 G gatkはJavaプログラムで、ローカルにダウンロードして解凍すれば使用できます.win 10でIDMを使用してgatk 4をダウンロードする.0.10.1アドレス格納ディレクトリ
/home/chaim/disk/gatk/
2.GATK4.0.10.1概要よく使われるpipelineは5種類あります Germline SNPs+Indels種系SNP+Indel Somatic SNVs+Indels体細胞単塩基突然変異 RNAseq SNPs + Indels Germline CNVs種系コピー数変異(Copy numbervariations,CNVs)は主に1 kb以上のDNA断片の欠失、挿入、重複などを指す.一般に構造的変異 Somatic CNVs体細胞コピー数変異1、2、4、5はDNAシーケンシング分析に適し、3はRNAシーケンシング分析に適している.公式文書 分析開始 GATK4.0全ゲノムと外現子群分析実戦
ソフトウェア: fastqc検出品質 fastq/trimmomatic品質管理 bwa対 samtoolフォーマット変換 データ格納場所
すべてのデータ環境の前提は、
1.品質管理検査
fastpのインストール
fastp品質管理データの使用
~~fastpはtrimmomaticより速く、効果が高いと言われています.一応信じる.
fastpパラメータ参照アドレス
bwaフローパラメータbwaリファレンス1 bwaリファレンス2
2.1. bwaインデックスファイルの作成
/bwaのコマンドはnohupを使用しないでください.Nohupの出力情報はbwaによってターゲットファイルに出力され、後続のステップ/
B 73シーケンスアドレス位置
/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa
/#注意:ここの2.2.1と2.2.2の2つのステップのbwaは必ず同じバージョンのbwaを使用しなければなりません.そうしないと、後でエラーが発生します.
2.2と2.3のどちらかを選択すればよい、2.3を使用することを提案する.bwaのmemの効率はより高く,より正確である.
2.2 bwa入力readsファイルのSA座標を探す
2.2.2 sai回転sam
注記:
2.3 BWAのmemの使用は、高速で一歩前進することができます.リファレンスアドレス(注意2.2と2.3、2つ選択すればよい、2.3を推奨)
3.Samファイルの並べ替え(recorder)
ダウンロードインストールの最新版picardは、パス
4.samファイルをbamファイルに変換します.
5.bamファイルのsortソート
6. Merge
7. Duplicates Marking
一般にRNA−seqはこのステップを必要とせず、chip−seq、dap−seq、call snpはこのステップを必要とする.具体的な原理シーケンシング原理はランダムに遮断され,理論的に2つの完全に同じreadが現れる確率は非常に低く,ライブラリの構築時にPCR増幅に偏向性があるため,完全に同じreadを示す.
8.前のステップの結果を生成するインデックスファイル
/#前のbwaのmemのRパラメータのため、私は第1回の実行時に完全に設定していないので、ここでヘッドファイルを2回変更する必要があります*/picardを使用してヘッドファイルID strを変更します:readsセットID番号を入力します;LB:readセットライブラリ名;PL:シーケンシングプラットフォーム(illunimaまたはsolid);PU:シーケンシングプラットフォームの下位単位名(runの名称);SM:サンプル名.
/必ず手動でヘッダーファイルを追加しないでください.手動の後続は認識できません./
9.Base (Quality Score) Recalibration
Tools involved:BaseRecalibrator,Apply Recalibration,AnalyzeCovariates(optional)参照アドレスフロー参照アドレス塩基質量分数再較正(Base quality score recalibration,BQSR)は,機械学習により元の塩基の質量分数を調整する.2つのステップに分けられます.既存のsnpデータベースを用いて相関モデルを構築し、再較正テーブル を生成する.は、このモデルに従って元の塩基を調整し、非既知のSNP領域のみを調整する.パラメータリスト-R:参照ゲノム-I:入力BAMファイル--known-sites既知SNPのvcfファイル-O:出力リキャリブレーションテーブル
#上記で生成したbamファイルが使用可能かどうかを検出します.
二、GATK変異検出
チュートリアルのアドレスの説明を参照してください.その後、24プロセスを使用することを表す
1.2以前に生成されたGVCFファイルを1つのファイルにマージする.
1.3 GenotypeGVCFs
2.select SNP
3.select indel
/#4.1および4.2変異フィルタリングは、異なるアルゴリズムのフィルタリングである.4.1は機械パラメータフィルタリング、4.2は機械学習フィルタリングである.どちらかを選択すると/4.1 filter SNP(変異フィルタ、ハードフィルタ)になります.パラメータの説明
パラメータの詳細(参照)
以下のフィルタパラメータは,自分の実際のvcfの値に基づいて分布し,自分で修正することができる.
QDは,単位深さの変異値を記述し,大きいほど信頼性が高い.一般に<2の値をフィルタします.FSは、正負鎖の特異性を記述し、差異が大きく、シーケンシングまたは組立の過程でランダムではないことを示している.FSは小さいほど良い.一般に、bwa−memを使用する場合、正常値は60であるべきであり、ある部位シーケンシングreadsの質量値の離散度を記述する.MQ<40.0 MQRankSum<-12.5 SORも、正負鎖特異性を表し、正常値が0-3、フィルタリング>3の値である.
4.2変種質量制御VQSRは、2つのステップに分けられる(##このステップは本実験では適用されず、実行されていない.)
/本実験では,このモデルを用いてフィルタリングすることはできず,このモデルは多様なvcf質量制御に適応している/ VariantRecalibratorは、フィルタリング(VariantRecalibrator) のためにバリアント品質を評価するための再較正モデルを構築する.変異質量スコア再較正ApplyVQSR
VariantRecalibrator
VQSR変換vcfはtable形式 である.
三、qtl-seqrパッケージ(第4ステップと同時に行うことができる)QTLseqrをgithubで使用する
以下のコードはR言語コード
四.gatk 4出力のvcfはsnpEffを用いて処理される.snpEffのチュートリアルダウンロードsnpEffアドレス解凍 ですトウモロコシzm 437バージョンを構成するデータベースは、パス を作成します.
上記のゲノム注釈ファイルはgff 3フォーマットであることに注意してください.snpEffを修正します.configのパラメータは次のように追加されます.
snpEFFディレクトリに戻り、コマンド vcf形式のファイルについて注記: が実行する.
3つのファイルが出力されますsnpEff_genes.txt snpEff_summary.html common.eff.vcf他の注釈ファイルの使用を変更したい場合は、
GoogleアカウントGATK 3のダウンロードアドレスGATK 4の各バージョンにログインする必要があり、コマンドによっては変更があります.1.GATKのインストール、他人のチュートリアルテンセントクラウドを使用するチュートリアルGATK 4データ前処理変異検出(BWA+SAMtools+picard+GATK)bwa+samtools+picardtools+GATK call SNPプロセスは、私たちのサーバーのスラグのネットワーク問題のため、内容はすべて地元のwin 10にダウンロードした後、サーバーにアップロードします.サーバ構成:cpu 4個6コア2スレッド=48メモリ:396 G gatkはJavaプログラムで、ローカルにダウンロードして解凍すれば使用できます.win 10でIDMを使用してgatk 4をダウンロードする.0.10.1アドレス格納ディレクトリ
/home/chaim/disk/gatk/
unzip gatk4.0.10.1
cd gatk4.0.10.1/
chmod 777 gatk
./gatk --list // gatk
2.GATK4.0.10.1概要よく使われるpipelineは5種類あります
ソフトウェア:
すべてのデータ環境の前提は、
/home/chaim/disk/BSA/
ディレクトリにディレクトリファイルがあります.119-8-1 //119-8 1
-A23-16551278-1279119-8_combined_R1.fastq.gz
-A23-16551278-1279119-8_combined_R2.fastq.gz
119-8-2 //119-8 2
-A23-16551278-1279-119-8_combined_R1.fastq.gz
-A23-16551278-1279-119-8_combined_R2.fastq.gz
origin
- B17SF2447-20_L1_358358.R1.clean.fastq_2.gz
- B17SF2447-20_L2_358358.R1.clean.fastq.gz
- B17SF2447-20_L1_358358.R2.clean.fastq.gz
- B17SF2447-20_L2_358358.R2.clean.fastq.gz
//
1.品質管理検査
fastqc *.fastq.gz -t 8 -o fastqc_out/
fastpのインストール
wget http://opengene.org/fastp/fastp
chmod 755 fastp
fastp品質管理データの使用
~~fastpはtrimmomaticより速く、効果が高いと言われています.一応信じる.
./fastp -i in.R1.fq.gz -o out.R1.fq.gz -I in.R2.fq.gz -O out.R2.fq.gz
# /BSA/119-8 119-8
../origin/fastp -i ../119-8-1/A23-16551278-1279119-8_combined_R1.fastq.gz -o ./fastp_out/119-8_1.R1.clean.fastq.gz -I ../119-8-1/A23-16551278-1279119-8_combined_R2.fastq.gz -O ./fastp_out/119-8_1.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &
../origin/fastp -i ../119-8-2/A23-16551278-1279-119-8_combined_R1.fastq.gz -o ./fastp_out/119-8_2.R1.clean.fastq.gz -I ../119-8-2/A23-16551278-1279-119-8_combined_R2.fastq.gz -O ./fastp_out/119-8_2.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &
# /BSA/origin/
./fastp -i ./B17SF2447-20_L1_358358.R1.clean.fastq.gz -o ./fastp_out/2447-20_L1$.R1.clean.fastq.gz -I ./B17SF2447-20_L1_358358.R2.clean.fastq.gz -O ./fastp_out/2447-20_L1.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &
./fastp -i ./B17SF2447-20_L2_358358.R1.clean.fastq.gz -o ./fastp_out/2447-20_L2.R1.clean.fastq.gz -I ./B17SF2447-20_L2_358358.R2.clean.fastq.gz -O ./fastp_out/2447-20_L2.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &
fastpパラメータ参照アドレス
-i R1
両端シーケンシングデータを入力R 1端-o outputR1
質量制御後出力R 1端-I R2
入力R 2原始シーケンシングデータ-O outputR2
質量制御後出力R 2端-Q
質量フィルタ--thread=5
無効設定スレッド数5 --length_required=50
設定フィルタの最短シーケンス長50 bp --n_base_limit=6
1つのreadsにおけるNの回数が6より大きく、このreads --compression=6
出力を破棄したgzipファイルの圧縮度合いは6,1−9であり、圧縮度合いは増大する.bwaフローパラメータbwaリファレンス1 bwaリファレンス2
2.1. bwaインデックスファイルの作成
/bwaのコマンドはnohupを使用しないでください.Nohupの出力情報はbwaによってターゲットファイルに出力され、後続のステップ/
B 73シーケンスアドレス位置
/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa
/BSA/bwa/zm437
// /BSA/bwa/
bwa index -a bwtsw -p zm437 /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa
index -a bwtsw
設定モード、大ゲノム-p zm437
設定出力ファイル名に適合/#注意:ここの2.2.1と2.2.2の2つのステップのbwaは必ず同じバージョンのbwaを使用しなければなりません.そうしないと、後でエラーが発生します.
2.2と2.3のどちらかを選択すればよい、2.3を使用することを提案する.bwaのmemの効率はより高く,より正確である.
2.2 bwa入力readsファイルのSA座標を探す
// /BSA/bwa/
bwa aln zm437 read1.fq.gz -l 30 -k 2 -t 8 -I > read1.fq.gz.sai
bwa aln zm437 read2.fq.gz -l 30 -k 2 -t 8 -I > read2.fq.gz.sai
// 4 2447-20
bwa aln zm437 ../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L1.R1.fq.gz.sai &
bwa aln zm437 ../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L1.R2.fq.gz.sai &
bwa aln zm437 ../origin/fastp_out/2447-20_L2.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L2.R1.fq.gz.sai &
bwa aln zm437 ../origin/fastp_out/2447-20_L2.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L2.R2.fq.gz.sai &
// 4 119-8
bwa aln zm437 ../119-8/fastp_out/119-8_1.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_1.R1.fq.gz.sai &
bwa aln zm437 ../119-8/fastp_out/119-8_1.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_1.R2.fq.gz.sai &
bwa aln zm437 ../119-8/fastp_out/119-8_2.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_2.R1.fq.gz.sai &
bwa aln zm437 ../119-8/fastp_out/119-8_2.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_2.R2.fq.gz.sai &
2.2.2 sai回転sam
bwa sampe -r "@RG\tID:\tLB:\tSM:\tPL:ILLUMINA" read1.fq.gz.sai read2.fq.gz.sai read1.fq.gz read2.fq.gz > read.sam
注記:
SAMPLE_NAME
は、対応するサンプル名に置き換える必要があります.そうしないと、サンプルとして処理されます.//2447-20
bwa sampe zm437 -r "@RG\tID:2447-20\tLB:B73\tSM:2447-20_L1\tPL:ILLUMINA" 2447-20_L1.R1.fq.gz.sai 2447-20_L1.R2.fq.gz.sai ../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz >2447-20_L1.sam &
bwa sampe zm437 -r "@RG\tID:2447-20\tLB:B73\tSM:2447-20_L2\tPL:ILLUMINA" 2447-20_L2.R1.fq.gz.sai 2447-20_L2.R2.fq.gz.sai ../origin/fastp_out/2447-20_L2.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L2.R2.clean.fastq.gz >2447-20_L2.sam &
//119-8
bwa sampe zm437 -r "@RG\tID:119-8\tLB:B73\tSM:119-8_1\tPL:ILLUMINA" 119-8_1.R1.fq.gz.sai 119-8_1.R2.fq.gz.sai ../119-8/fastp_out/119-8_1.R1.clean.fastq.gz ../119-8/fastp_out/119-8_1.R2.clean.fastq.gz >119-8_1.sam &
bwa sampe zm437 -r "@RG\tID:119-8\tLB:B73\tSM:119-8_2\tPL:ILLUMINA" 119-8_2.R1.fq.gz.sai 119-8_2.R2.fq.gz.sai ../119-8/fastp_out/119-8_2.R1.clean.fastq.gz ../119-8/fastp_out/119-8_2.R2.clean.fastq.gz >119-8_2.sam &
2.3 BWAのmemの使用は、高速で一歩前進することができます.リファレンスアドレス(注意2.2と2.3、2つ選択すればよい、2.3を推奨)
#bwa mem
/* /home/chaim/disk/BSA/bwa/*/
/*zm437 B73 */
/* -R */
bwa mem -t 24 -M -P -R '@RG\tID:2447-20\tSM:2447-20\tLB:WES\tPL:Illumina' zm437 ../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz >2447-20_L1.sam &
bwa mem -t 24 -M -P -R '@RG\tID:2447-20\tSM:2447-20\tLB:WES\tPL:Illumina' zm437 ../origin/fastp_out/2447-20_L2.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L2.R2.clean.fastq.gz >2447-20_L2.sam &
bwa mem -t 12 -M -P -R '@RG\tID:119-8\tSM:119-8\tLB:WES\tPL:Illumina' zm437 ../119-8/fastp_out/119-8_1.R1.clean.fastq.gz ../119-8/fastp_out/119-8_1.R2.clean.fastq.gz >119-8_1.sam &
bwa mem -t 8 -M -P -R '@RG\tID:119-8\tSM:119-8\tLB:WES\tPL:Illumina' zm437 ../119-8/fastp_out/119-8_2.R1.clean.fastq.gz ../119-8/fastp_out/119-8_2.R2.clean.fastq.gz >119-8_2.sam &
3.Samファイルの並べ替え(recorder)
ダウンロードインストールの最新版picardは、パス
/home/chaim/disk/BSA/bwa/
に保存され、このパスでjava -jar picard.jar -h
が実行され、picardに含まれるすべてのツールがリストされます.3.1インデックスシーケンスの構築nohup samtools faidx zm437 &
3.2 Samファイルの並べ替えjava -jar picard.jar CreateSequenceDictionary R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa O=zm437.dict
java -jar picard.jar ReorderSam I=2447-20_L1.sam O=2447-20_L1.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa &
java -jar picard.jar ReorderSam I=2447-20_L2.sam O=2447-20_L2.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa &
java -jar picard.jar ReorderSam I=119-8_1.sam O=119-8_1.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa
java -jar picard.jar ReorderSam I=119-8_2.sam O=119-8_2.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa &
4.samファイルをbamファイルに変換します.
samtools view --threads 24 -bS 2447-20_L1.reordered.sam -o 2447-20_L1.bam &
samtools view --threads 24 -bS 2447-20_L2.reordered.sam -o 2447-20_L2.bam &
samtools view --threads 8 -bS 119-8_1.reordered.sam -o 119-8_1.bam
samtools view --threads 8 -bS 119-8_2.reordered.sam -o 119-8_2.bam &
5.bamファイルのsortソート
java -jar picard.jar SortSam INPUT=2447-20_L1.bam OUTPUT=2447-20_L1.sort.bam SORT_ORDER=coordinate &
java -Xmx48G -jar picard.jar SortSam INPUT=2447-20_L2.bam OUTPUT=2447-20_L2.sort.bam SORT_ORDER=coordinate &
java -Xmx96G -jar picard.jar SortSam INPUT=119-8_1.bam OUTPUT=119-8_1.sort.bam SORT_ORDER=coordinate
java -jar picard.jar SortSam INPUT=119-8_2.bam OUTPUT=119-8_2.sort.bam SORT_ORDER=coordinate &
6. Merge
\\ lane bam 。
java -jar picard.jar MergeSamFiles I=2447-20_L1.sort.bam I=2447-20_L2.sort.bam O=2447-20.bam &
java -jar picard.jar MergeSamFiles I=119-8_1.sort.bam I=119-8_2.sort.bam O=119-8.bam
7. Duplicates Marking
一般にRNA−seqはこのステップを必要とせず、chip−seq、dap−seq、call snpはこのステップを必要とする.具体的な原理シーケンシング原理はランダムに遮断され,理論的に2つの完全に同じreadが現れる確率は非常に低く,ライブラリの構築時にPCR増幅に偏向性があるため,完全に同じreadを示す.
java -jar picard.jar MarkDuplicates REMOVE_DUPLICATES= false MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 INPUT=2447-20.bam OUTPUT=2447-20.repeatmark.bam METRICS_FILE=2447-20.bam.metrics
java -jar picard.jar MarkDuplicates REMOVE_DUPLICATES= false MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 INPUT=119-8.bam OUTPUT=119-8.repeatmark.bam METRICS_FILE=119-8.bam.metrics
8.前のステップの結果を生成するインデックスファイル
samtools index 2447-20.repeatmark.bam
samtools index 119-8.repeatmark.bam
/#前のbwaのmemのRパラメータのため、私は第1回の実行時に完全に設定していないので、ここでヘッドファイルを2回変更する必要があります*/picardを使用してヘッドファイルID strを変更します:readsセットID番号を入力します;LB:readセットライブラリ名;PL:シーケンシングプラットフォーム(illunimaまたはsolid);PU:シーケンシングプラットフォームの下位単位名(runの名称);SM:サンプル名.
java -Xmx96g -jar picard.jar AddOrReplaceReadGroups I=2447-20.repeatmark.bam O=2447-20.repeat.bam LB=lib2447-20 PL=illumina PU=2447-20 SM=2447-20 &
java -Xmx96g -jar picard.jar AddOrReplaceReadGroups I=119-8.repeatmark.bam O=119-8.repeat.bam LB=lib119-8 PL=illumina PU=119-8 SM=119-8 &
/必ず手動でヘッダーファイルを追加しないでください.手動の後続は認識できません./
9.Base (Quality Score) Recalibration
Tools involved:BaseRecalibrator,Apply Recalibration,AnalyzeCovariates(optional)参照アドレスフロー参照アドレス塩基質量分数再較正(Base quality score recalibration,BQSR)は,機械学習により元の塩基の質量分数を調整する.2つのステップに分けられます.
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar BaseRecalibrator -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 2447-20.repeat.bam --known-sites /home/guo/maize/zm437/zea_mays_vcfsort.vcf -O 2447-20_recal_data.table &
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyBQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 2447-20.repeat.bam -bqsr 2447-20_recal_data.table -O 2447-20_recal_reads.bam &
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar BaseRecalibrator -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 119-8.repeat.bam --known-sites /home/guo/maize/zm437/zea_mays_vcfsort.vcf -O 119-8_recal_data.table &
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyBQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 119-8.repeat.bam -bqsr 119-8_recal_data.table -O 119-8_recal_reads.bam &
#上記で生成したbamファイルが使用可能かどうかを検出します.
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ValidateSamFile -I 2447-20_recal_reads.bam
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ValidateSamFile -I 119-8_recal_reads.bam
no errors found
が表示する場合、HaplotyCaller call SNP/Indelを使用することができる.二、GATK変異検出
チュートリアルのアドレスの説明を参照してください.その後、24プロセスを使用することを表す
-nt 24
というパラメータが一部のコマンドで使用されます.各コマンドがマルチプロセスを開くわけではありません.gatk公式サイトでドキュメントを検索し、コマンドを検索した後、コマンドのAPIドキュメントでthread
を検索すると、マルチスレッドが使用できるかどうかをすばやく検索できます.(このステップV 4.1.4.0ではjavaのメモリ割り当ては16 gを超えないことが望ましいが、8 gで十分であることを提案する.hmmプロセスもポイント4-8個を少なく設定する.Javaメモリ設定が高すぎると、サーバのメモリ範囲でもエラーが発生し、メモリ割り当てが不足しているため途中でプログラムを削除し、途中で切ると同様に修復方法がある.picardは中断したgvcfファイルを修復する).raw vcfファイルを生成します(3-5日で実行できるのは、すべて大きなサーバです)パラメータは、HaplotypeCaller
でgvcfファイルを生成してからCombineGVCFs
を実行することを示します.java -Xmx96G -jar/home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.JAr#Xmx 96 Gで使用する最大メモリHaplotyCaller#HaplotyCallerモードを使用し、食事構成-R/home/chaim/disk/BSA/bwa/zm 437#参照B 73ゲノム-I 2447-20.repeatmark.Bam#複数サンプルの場合-I sample 1.bam -I sample2.bam--dbsnp zm 437 vcf#B 73参照snp-stand_emit_conf 10 -stand_call_conf 30 -O 2447-20.vcf 1.1 Gvcfファイルを生成する(このステップは時間の無駄で、3-5昼夜以上かかる)java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar HaplotypeCaller -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 119-8_recal_reads.bam --dbsnp /home/guo/maize/zm437/zea_mays_vcfsort.vcf --native-pair-hmm-threads 96 -stand-call-conf 30 -O 119-8.g.vcf.gz -ERC GVCF
java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar HaplotypeCaller -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 2447-20_recal_reads.bam --dbsnp /home/guo/maize/zm437/zea_mays_vcfsort.vcf --native-pair-hmm-threads 96 -stand-call-conf 30 -O 2447-20.g.vcf.gz -ERC GVCF
1.2以前に生成されたGVCFファイルを1つのファイルにマージする.
java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar CombineGVCFs -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa --dbsnp /home/guo/maize/zm437/zea_mays_vcfsort.vcf --variant 2447-20.g.vcf.gz --variant 119-8.g.vcf.gz -O cohort.g.vcf.gz
1.3 GenotypeGVCFs
java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar GenotypeGVCFs -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V cohort.g.vcf.gz -O common.vcf.gz
2.select SNP
java -Xmx96g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar SelectVariants -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -O common_SNP.vcf --variant common.vcf.gz --select-type-to-include SNP 2>select_snp.err
3.select indel
java -Xmx96g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar SelectVariants -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -O common_INDEL.vcf --variant common.vcf.gz --select-type-to-include INDEL 2>select_indel.err
/#4.1および4.2変異フィルタリングは、異なるアルゴリズムのフィルタリングである.4.1は機械パラメータフィルタリング、4.2は機械学習フィルタリングである.どちらかを選択すると/4.1 filter SNP(変異フィルタ、ハードフィルタ)になります.パラメータの説明
java -Xmx4g -jar $GATK -R $REF -T VariantFiltration --variant $Slect_SNP --clusterSize 4 --clusterWindowSize 10 --maskName aroundIndel --mask $Slest_INdel -maskExtend 3 --filterName "lowMQRankSum" --filterExpression "MQRankSum < -12.5" --filterName "highFS" --filterExpression "FS > 60.0" --filterName "lowReadPosRankSum" --filterExpression "ReadPosRankSum < -8.0" --filterName "lowMQ" --filterExpression "MQ < 40.0" --filterName "lowQD" --filterExpression "QD < 2.0" --out $Filter_SNP --genotypeFilterName "lowDP" --genotypeFilterExpression "DP < 8.0" >filter_snp.err
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar VariantFiltration -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa --variant common_SNP.vcf --cluster-size 4 --cluster-window-size 10 --mask-name aroundIndel --mask common_INDEL.vcf -mask-extension 3 --filter-name "lowMQRankSum" --filter-expression "QUAL < 30" --filter-name "qua130" --filter-expression "MQRankSum < -12.5" --filter-name "highFS" --filter-expression "FS > 60.0" --filter-name "lowReadPosRankSum" --filter-expression "ReadPosRankSum < -8.0" --filter-name "lowMQ" --filter-expression "MQ < 40.0" --filter-name "lowQD" --filter-expression "QD < 2.0" -O common_filtration.vcf --genotype-filter-name "lowDP" --genotype-filter-expression "DP < 8.0" >filter_snp.err
パラメータの詳細(参照)
以下のフィルタパラメータは,自分の実際のvcfの値に基づいて分布し,自分で修正することができる.
QDは,単位深さの変異値を記述し,大きいほど信頼性が高い.一般に<2の値をフィルタします.FSは、正負鎖の特異性を記述し、差異が大きく、シーケンシングまたは組立の過程でランダムではないことを示している.FSは小さいほど良い.一般に、bwa−memを使用する場合、正常値は60であるべきであり、ある部位シーケンシングreadsの質量値の離散度を記述する.MQ<40.0 MQRankSum<-12.5 SORも、正負鎖特異性を表し、正常値が0-3、フィルタリング>3の値である.
4.2変種質量制御VQSRは、2つのステップに分けられる(##このステップは本実験では適用されず、実行されていない.)
/本実験では,このモデルを用いてフィルタリングすることはできず,このモデルは多様なvcf質量制御に適応している/
VariantRecalibrator
gatk VariantRecalibrator \
-R Homo_sapiens_assembly38.fasta \
-V input.vcf.gz \
--resource hapmap,known = false,training = true,truth = true,prior = 15.0:hapmap_3.3.hg38.sites.vcf.gz \
--resource omni,known = false,training = true,truth = false,prior = 12.0:1000G_omni2.5.hg38.sites.vcf.gz \
--resource 1000G,known = false,training = true,truth = false,prior = 10.0:1000G_phase1.snps.high_confidence.hg38.vcf.gz \
--resource dbsnp,known = true,training = false,truth = false,prior = 2.0:Homo_sapiens_assembly38.dbsnp138.vcf.gz \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
-mode SNP
-O output.recal \
--tranches-file output.tranches \
--rscript-file output.plots.R
VQSR
gatk ApplyVQSR \
-R Homo_sapiens_assembly38.fasta \
-V input.vcf.gz \
-O output.vcf.gz \
--truth-sensitivity-filter-level 99.0 \
--tranches-file output.tranches \
--recal-file output.recal \
-mode SNP
java -Xmx128g jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyVQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V 2447-20.vcf -O 2447-20.vqsr.vcf &
java -Xmx128g jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyVQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V 119-8.vcf -O 119-8.vqsr.vcf &
java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar VariantsToTable -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V common_filtration.vcf -F CHROM -F POS -F REF -F ALT -GF AD -GF DP -GF GQ -GF PL -O common.table
三、qtl-seqrパッケージ(第4ステップと同時に行うことができる)QTLseqrをgithubで使用する
以下のコードはR言語コード
# QTL-seqr
#
installed.packages('devtools')
library('devtools')
install_github("bmansfeld/QTLseqr")
#load the package
library("QTLseqr")
#set sample and file name
HighBulk
四.gatk 4出力のvcfはsnpEffを用いて処理される.snpEffのチュートリアル
unzip snpEff_latest_core.zip
私の経路は/home/chaim/bsa/snpEff//home/chaim/bsa/snpEff/snpEff/
ディレクトリの下にフォルダdata
、cd /home/chaim/bsa/snpEff/snpEff/
mkdir data
cd data
mkdir genomes
mkdir zm437
# zm437 genes.gff, ,protein.fa
# genomes zm437.fa
上記のゲノム注釈ファイルはgff 3フォーマットであることに注意してください.snpEffを修正します.configのパラメータは次のように追加されます.
#maize genome,version zm437
zm437.genome:maize
snpEFFディレクトリに戻り、コマンド
java -jar snpEff.jar build -gff3 -v zm437
を実行bwa
ディレクトリにはGATK 4処理後のファイルcommon_filtration.vcf
が格納されている/home/chaim/bsa/bwa
ディレクトリでは、次のコマンド java -Xmx128g -jar /home/chaim/bsa/snpEff/snpEff/snpEff.jar zm437 common_filtration.vcf > common.eff.vcf
3つのファイルが出力されますsnpEff_genes.txt snpEff_summary.html common.eff.vcf
/home/chaim/bsa/snpEff/snpEff//data/zm437/snpEffectPredictor.bin
を削除して削除すればよい.手順1からやり直してください.