全基因组测序流程 | 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文件

1	`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

进入R里面合并文件

参考：http://localhost:17435/notebooks/projects/TAP_WGS/analysis/RNA-seq.ipynb

从gtf文件里提取信息

grep -v '#' Homo_sapiens_assembly19.gtf | cut -f9 | cut -d' ' -f2,10 | uniq > id.symbol
 
grep -v '#' Homo_sapiens_assembly19.gtf | cut -f1-5 > id.length
grep -v '#' Homo_sapiens_assembly19.gtf | cut -f9 | cut -d' ' -f2,10 > full.id.symbol

转化count为TPM

# counts:转录组的count矩阵，行为基因，列为样本
# effLen：一个数值型向量，值是基因长度，顺序应该与count的列一致对应。
Counts2TPM <- function(counts, effLen){
  rate <- log(counts) - log(effLen)
  denom <- log(sum(exp(rate)))
  exp(rate - denom + log(1e6))
}