全基因组测序流程 | WGS pipeline
创建conda环境，安装必要软件
conda create -n wgs
conda activate wgs
conda install bioconda::bwa
下载最佳ref fasta
gcloud storage cp gs://BUCKET_NAME/OBJECT_NAME SAVE_TO_LOCATION
gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.fasta
Resource bundle
https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0/
构建索引【可以直接下载，broad有标准流程】
../refdata-gex-GRCh38-2020-A/fasta/genome.fa
bwa index -p GRCh38 hg38.fasta
/home/zz950/reference/bwa_wgs/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta.64
bwa标准比对
然后合并两个lane
部分下载的数据以CRAM格式保存
How to convert CRAM back to BAM using samtools
It's the same command, just use the -b option instead of -C.
需要最新版samtools读取cram文件
conda install bioconda::samtools
下载GATk： https://github.com/broadinstitute/gatk/releases
下载Jar：https://gatk.broadinstitute.org/hc/en-us/community/posts/360073027611-Where-and-how-I-can-get-GenomeAnalysisTK-jar-script
VariantRecalibrator java.lang.IllegalArgumentException: No data found.
No data found when using VariantRecalibrator
Using CRAM within Samtools
还有配对的RNA-seq
但是使用的是hg19的参考基因组，用的STAR比对。
https://data.broadinstitute.org/snowman/hg19/annotation/
这个bam可以用htseq-count来计数
conda install -c bioconda htseq
conda install bioconda::htseq
##合并htdeq生成的count文件到matrix.count
cd $MATRIX/htseqperl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count | grep -v "_" >matrix.count
##合并featureCounts生成的count文件到matrix_2.count
cd $MATRIX/featureCountsfor i in `seq 15 18`;docut -f 1,7 SRR20382${i}.feature.count | grep -v "^#" > SRR20382${i}.matrixsed -i '1d' SRR20382${i}.matrixdoneperl -lne 'if ($ARGV=~/(.*).matrix/){print "$1\t$_"}' *.matrix >matrix_2.count
参考流程：
- 从零开始完整学习全基因组测序（WGS）数据分析之构建WGS主流程
- 使用htseq-count进行定量分析