2021/03/25
一、測序數(shù)據(jù)下載
NCBI SRA 數(shù)據(jù)庫
sra
wget ftp://ftp.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-centos_linux64.tar.gz
tar -zxvf sratoolkit.current-centos_linux64.tar.gz
wget https://ftp.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
tar -zxvf sratoolkit.current-ubuntu64.tar.gz
./prefetch SRR13948848 # 下載數(shù)據(jù)
md5sum -b # md5校驗(yàn)蝶怔，驗(yàn)證數(shù)據(jù)有沒有下載完全
./fastq-dump -X 5 -Z SRR13948848 # 將sra轉(zhuǎn)換為fastq
【使用說明】
二、fastqc （fastq文件）
（一）瀑凝、下載和安裝
下載
wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.9.zip
unzip fastqc_v0.11.9.zip
chmod 755 fastqc
export PATH=/home/user/FastQC/:$PATH
source ~/.bashrc
fastqc --help
（二）、在 python 中 使用
先寫出 需要的 命令行
-o --outdir 輸出路徑
-f --format 指定輸入文件的格式
--extract 生成的報(bào)告默認(rèn)會(huì)打包成1個(gè)壓縮文件，使用這個(gè)參數(shù)是讓程序不打包
-t --threads 選擇程序運(yùn)行的線程數(shù)瞬场，每個(gè)線程會(huì)占用250MB內(nèi)存
--min_length 設(shè)置序列的最小長度
-c --contaminants 污染物選項(xiàng)
Sequence 可能的污染序列
-a --adapters adpater序列信息乾巧，默認(rèn)通用引物
-q --quiet 安靜運(yùn)行
import os
commond = f'{fastqc_path} --extract -q -o {output_dir} {fastq1_path} {fastq2_path}'
os.system(commond)
生成文件：html 和 zip 文件
三句喜、blat （fasta文件 ------> psl 文件）
下載
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64.v385/blat/blat
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/blat/blat
用法
blat database query [-ooc=11.ooc] output.psl
輸入文件.fa, .nib, .2bit file
-t=type 數(shù)據(jù)庫類型：dna(default) / prot(protein) / dnax(6種）
-q=type 輸入類型：dna / rna / prot / dnax / rnax
-prot Synonymous with -t=prot -q=prot
-minIdentity=N 最小序列相似度，默認(rèn)dna/rna 90,蛋白沟于，轉(zhuǎn)錄本25
-out=type 輸出文件格式：psl(默認(rèn)) / pslx / axt / maf / sim4 /
wublast / blast / blast8 / blast9
commond = f'{blat_software} -minIdentity=80 -out=psl {database} {fasta_file} {blat_psl}'
針對psl文件表現(xiàn)出的比對情況咳胃，進(jìn)行整理，過濾（同源性低reads旷太，剪切掉引物等）
''' blat_psl file
match mismatch repmatch N's Q_gap_count Q_gap_bases T_gap_count T_gap_bases
strand Q_name Q_size Q_start Q_end T_name T_size T_start T_end
block_count blocksizes qStarts tStarts
'''
四展懈、bwa + sambamba （fastq文件 ------> bam 文件）
wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.17.tar.bz2
tar jxvf bwa-0.7.17.tar.bz2
cd bwa-0.7.17
make
BWA是一個(gè)軟件包，用于將低差異序列映射到大型參考基因組（例如人類基因組）上供璧，由三種算法組成：BWA-backtrack, BWA-SW 和 BWA-MEM存崖。
bwa mem [ -aCHMpP ] [ -t nThreads ] [ -k minSeedLen ] [ -w bandWidth ] [ -d zDropoff ]
[ -r seedSplitRatio ] [ -c maxOcc ] [ -A matchScore ] [ -B mmPenalty ] [ -O gapOpenPen ]
[ -E gapExtPen ] [ -L clipPen ] [ -U unpairPen] [ -R RGline ] [ -v verboseLevel ]
db.prefix reads.fq [ mates.fq ]
-t 線程數(shù)[1] *
-k 最小種子長度[19]
-w 帶寬[100]
-d 離對角線的X-dropoff（Z-dropoff)[100]
-r 值越大，產(chǎn)生的種子越少睡毒，導(dǎo)致對齊速度更快来惧，但準(zhǔn)確性較低。[1.5]
-c 如果MEM在基因組中發(fā)生的次數(shù)超過INT演顾，則將其丟棄 供搀。[10000]
-P 在配對端模式下隅居，執(zhí)行SW僅可挽救丟失的匹配，但不要嘗試查找適合合適的匹配的匹配葛虐。
-A 匹配分?jǐn)?shù)胎源。[1]
-B 不匹配的罰款[4]
-O 差距開放處罰[6]
-E 間隙延長罰分[1]
-L 罰款[5]
-U 未配對讀對的處罰。[9]
-p 假定第一個(gè)輸入查詢文件是交錯(cuò)的成對配對的FASTA / Q屿脐。
-R 完成讀取組標(biāo)題行涕蚤。'\t'可以在STR中使用， 并將在輸出SAM中轉(zhuǎn)換為TAB摄悯。讀取組ID將附加到輸出中的每個(gè)讀取赞季。
一個(gè)示例是"@RG\tID:foo\tSM:bar"。[空值] *
-T 不要輸出分?jǐn)?shù)低于INT的對齊方式 [30] *
-a 輸出所有發(fā)現(xiàn)的單端或未配對的成對末端讀段的比對奢驯。這些比對將被標(biāo)記為次要比對申钩。
-C 將FASTA / Q注釋附加到SAM輸出。此選項(xiàng)可用于將讀取的元信息（例如條形碼）傳輸?shù)絊AM輸出瘪阁。
請注意撒遣，F(xiàn)ASTA / Q注釋（標(biāo)題行中的空格之后的字符串）必須符合SAM規(guī)范（例如BC：Z：CGTAC）。
格式錯(cuò)誤的注釋會(huì)導(dǎo)致不正確的SAM輸出管跺。
-H 在SAM輸出中使用硬限幅“ H”义黎。映射長重疊群或BAC序列時(shí)，此選項(xiàng)可能會(huì)大大減少輸出的冗余豁跑。
-M 將較短的拆分匹配標(biāo)記為次要（以實(shí)現(xiàn)Picard兼容性）廉涕。 *
-v 控制輸出的詳細(xì)級別。整個(gè)BWA均未完全支持此選項(xiàng)艇拍。[3]
sambamba
git clone --recursive https://github.com/lomereiter/sambamba.git 好似不行
wget -c https://github.com/biod/sambamba/archive/v0.6.9.tar.gz
cd sambamba && make sambamba-ldmd2-64
【sambamba-view】 從SAM / BAM文件中提取信息的工具
sambamba view [OPTIONS] <input.bam|input.sam> [region1 [...]]
-F, --filter=FILTER 設(shè)置自定義過濾器
-f, --format=FORMAT 輸出格式 *
-h, --with-header 打印sam標(biāo)頭
-H狐蜕， --header 僅打印sam標(biāo)頭
-I， --reference-info 以JSON格式輸出stdout參考序列名稱和長度 *
-c, --count 僅將匹配記錄的數(shù)量輸出到stdout卸夕，-hHI選項(xiàng)將被忽略
-v, --valid 僅輸出有效reads
-S, --sam-input 指定輸入SAM文件层释，默認(rèn)情況下，該工具將BAM文件作為輸入 *
-p, --show-progress 在終端中顯示進(jìn)度條快集，僅適用于BAM文件贡羔，沒有指定區(qū)域，僅在讀取完整文件時(shí)有效
-l, --compression-level=COMPRESSION_LEVEL 設(shè)置BAM輸出的壓縮級別个初，范圍為0到9乖寒。
-o, --output-filename=FILENAME 指定輸出文件名（默認(rèn)情況下，所有內(nèi)容均寫入stdout）
-t, --nthreads=NTHREADS 線程數(shù)
【sambamba-sort】 用于對BAM文件進(jìn)行排序的工具
sambamba sort OPTIONS <input.bam>
-m, --memory-limit=LIMIT 設(shè)置已用內(nèi)存的上限 *
--tmpdir=/tmp 使用TMPDIR輸出排序的塊院溺。默認(rèn)行為是使用系統(tǒng)臨時(shí)目錄 *
-o, --out=OUTPUTFILE 輸出文件名宵统。如未提供，則結(jié)果寫入.sorted.bam文件中 *
-n, --sort-by-name 按讀取名稱排序，而不是進(jìn)行坐標(biāo)排序
-l, --compression-level=COMPRESSION_LEVEL 用于排序的BAM的壓縮級別马澈，從0到9
-u, --uncompressed-chunks 將排序的塊寫為未壓縮的BAM。默認(rèn)使用壓縮級別1
-p, --show-progress 顯示類似wget的進(jìn)度條（實(shí)際上弄息，其中兩個(gè)接一個(gè)痊班，第一個(gè)用于排序，然后另一個(gè)用于合并）
-t, --nthreads=NTHREADS 線程數(shù) *
注：在Linux系統(tǒng)中摹量，標(biāo)準(zhǔn)輸入/dev/stdin , 標(biāo)準(zhǔn)輸出/dev/stdout
sambamba 把 sam 文件可以輸入到 /dev/stdin 文件中涤伐，再將文件中的內(nèi)容整理輸出到 bam 文件中
bwa 可以和 sambamba 聯(lián)合起來用，將 bwa 的結(jié)果 sam文件直接用sambamba處理成 bam
bwa_commond = bwa mem -R "@RG\tID:{sample}\tSM{sample}\tLB:WGS\tPL:Illumina" ' genome_database fastq_file | sambamba view -S -f bam /dev/stdin |sambamba sort --tmpdir=/tmp -o bam_file /dev/stdin'
看情況添加其他選項(xiàng)缨称，包括線程數(shù)凝果，占用內(nèi)存啥的
五、samtools + bcftools （bam文件 ------> pileup/vcf 文件）
samtools
wget https://sourceforge.net/projects/samtools/files/samtools/1.12/samtools-1.12.tar.bz2/download
【samtools】 序列比對/圖譜（SAM）格式的實(shí)用程序
samtools view [*options*] *in.sam*|*in.bam*|*in.cram* [*region*...]
-b 以BAM格式輸出
-o 輸出到文件[stdout]
-u 輸出未壓縮的BAM,可節(jié)省壓縮/解壓縮時(shí)間睦尽，當(dāng)將輸出通過管道傳遞到另一個(gè)samtools命令時(shí)器净，該選項(xiàng)是首選的。
【samtools mpileup】 從對齊方式生成“ pileup”文本格式
samtools mpileup [-EB] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
-A当凡，-count-orphans 不要在變體調(diào)用中跳過異常的讀取對
-b山害，-bam-list 輸入的BAM文件列表，每行一個(gè)文件[null]
-B沿量，--no-BAQ 禁用基本比對質(zhì)量（BAQ）計(jì)算浪慌，BAQ是讀取堿基未對齊的Phred標(biāo)度概率
-Q，--min-BQ 最低基準(zhǔn)質(zhì)量
-q朴则，--min-MQ 最低mapping質(zhì)量
-f权纤，-fasta-ref FAID格式 的 faidx索引參考文件
...
Samtools mpileup仍然可以產(chǎn)生VCF和BCF輸出（帶有 -g或-u），但是不建議使用此功能乌妒，并且在將來的版本中將刪除該功能汹想。請使用 bcftools mpileup（2021年3月17日發(fā)行的samtools-1.12手冊頁）
例：call SNPs and short indels
samtools mpileup -uf ref.fa aln.bam | bcftools call -mv > var.raw.vcf
bcftools filter -s LowQual -e '%QUAL<20 || DP>100' var.raw.vcf > var.flt.vcf
獲取測序深度
samtools depth -r chr12:126073855-126073965 Ip.sorted.bam
chr12 126073855 5
chr12 126073856 15
bcftools
git clone --recurse-submodules git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
cd bcftools
# The following is optional:
# autoheader && autoconf && ./configure --enable-libgsl --enable-perl-filters
make
為了使用BCFtools插件，必須設(shè)置此環(huán)境變量并指向正確的位置：
export BCFTOOLS_PLUGINS=/path/to/bcftools/plugins
【bcftools】 用于變體調(diào)用和操縱VCF和BCF的實(shí)用程序
bcftools [--version|--version-only] [--help] [COMMAND] [OPTIONS]
六芥被、freebayes （bam文件 ------> vcf 文件）
git clone --recursive git://github.com/ekg/freebayes.git
make
freebayes [option] bam.files > vcf.files
-b --bam 將文件添加到要分析的BAM文件集
-f --fasta-reference 分析所需參考序列
-t --targets 對bed文件中目標(biāo)進(jìn)行分析
-A --cnv-map 從bed文件中讀取cnv mapping
-v --vcf 輸出vcf格式
-C --min-alternate-count 最少 alt 堿基數(shù)欧宜，默認(rèn)2
-F --min-alternate-fraction 最少 alt 堿基比例，雙倍體默認(rèn) 0.2
--genotype-qualities 在輸出的vcf文件中有GQ
--min-mapping-quality 默認(rèn)1
--min-base-quality 默認(rèn)0
...
freebayes -h 查看文檔
freebayes --fasta-reference genome_db --bam example.bam --min-alternate-count 2 --min-alternate-fraction 0.2 --min-base-quality 20 --min-mapping-quality 1 > result.vcf
七拴魄、annovar （vcf文件 ------>avinput ------> 不同文件冗茸，關(guān)于注釋）
下載需提供郵箱
tar xvfz annovar.latest.tar.gz
annotate_variation.pl 主程序，下載數(shù)據(jù)庫匹中，注釋
coding_change.pl 推斷蛋白質(zhì)序列
convert2annovar.pl 將多種格式轉(zhuǎn)為.avinput
retrieve_seq_from_fasta.pl 自行建立其他物種的轉(zhuǎn)錄本
table_annovar.pl 注釋程序
variants_reduction.pl 靈活定制過濾注釋流程
humandb 人類注釋數(shù)據(jù)庫
下載數(shù)據(jù)庫https://annovar.openbioinformatics.org/en/latest/user-guide/download/
perl annotate_variation.pl -downdb avdblist -buildver hg38 -webfrom annovar path/to/humandb
-format vcf4 可用于將VCF文件轉(zhuǎn)換為ANNOVAR輸入格式
convert2annovar.pl -format vcf4 example.vcf > ex2.avinput
chr start end ref alt
# 基于基因的注釋 https://annovar.openbioinformatics.org/en/latest/user-guide/gene/
生成兩個(gè)輸出文件: ex1.variant_function （包含所有變體的注釋）
ex1.exonic_variant_function（包含由于外顯子變體導(dǎo)致的氨基酸變化）
（要更改輸出文件名夏漱，請使用--outfile參數(shù)）
# 基于區(qū)域的注釋 https://annovar.openbioinformatics.org/en/latest/user-guide/region/
annotate_variation.pl --regionanno
annotate_variation.pl -downdb phastConsElements46way humandb/
輸出保存在phastConsElements46way文件
# 基于過濾器的注釋 https://annovar.openbioinformatics.org/en/latest/user-guide/filter/
annotate_variation.pl --filter -dbtype database_name --buildver genome_version -out 輸出文件開頭
已知的變體將與等位基因頻率一起寫入_dropped文件
沒有匹配的數(shù)據(jù)庫條目的變量將被寫入_filtered文件
八、Intervar
InterVar是用于臨床意義變異解釋的python腳本
同作者關(guān)于Intervar
python Intervar.py -b genome_version -i example.bed --input_type=AVinput -d /InterVar_HumanDB -o output/dir'
九顶捷、分析 （vcf 及其各項(xiàng)注釋）
1挂绰、突變頻譜圖
2、突變熱圖
3、基因組疤痕評分
4葵蒂、突變負(fù)荷
5交播、驅(qū)動(dòng)突變
6、罕見突變與融合基因檢測
7践付、突變過程及其突變信號
8秦士、功能富集分析
（1）基因集富集分析
（2）路徑分析和可視化
（3）網(wǎng)絡(luò)分析
基本分析邏輯
