Charpter6：數(shù)據(jù)的格式Data Formats
常用數(shù)據(jù)庫(kù)
- NCBI
- EBI
- Uniprot
- Phytozome
- Ensemble plants
- TAIR
- PlantGDB
- PlantTFDB 植物轉(zhuǎn)錄因子數(shù)據(jù)庫(kù)
這些是我目前常用的一些數(shù)據(jù)庫(kù)店读，可以看出來(lái)我是做植物的。其實(shí)植物領(lǐng)域還有許多重要的數(shù)據(jù)庫(kù)唉匾，更不必說(shuō)動(dòng)物的了搁进，各個(gè)物種都能拿出來(lái)做個(gè)數(shù)據(jù)庫(kù)。數(shù)據(jù)這么多逞敷，如何更加高效的從數(shù)據(jù)庫(kù)中實(shí)現(xiàn)你想要的目的是門學(xué)問(wèn)雅宾。而自己目前對(duì)這些數(shù)據(jù)庫(kù)知道的也只是初步的一些東西。其實(shí)除了NCBI外很多重要的數(shù)據(jù)庫(kù)都有一些很實(shí)用的功能鸭丛，但是卻缺少教程指導(dǎo)推廣，以后有機(jī)會(huì)可以再對(duì)一些數(shù)據(jù)庫(kù)的功能進(jìn)行探索探索唐责。
數(shù)據(jù)格式
1. GeneBank format
- 最常見(jiàn)的生物數(shù)據(jù)存儲(chǔ)方式FASTA,GeneBank, FASTQ,前兩者被稱為curated sequence information.
- Genebank/gb格式文件為NCBI存儲(chǔ)基因信息的格式鳞溉。其中NCBI里的RefSeq數(shù)據(jù)庫(kù)為標(biāo)準(zhǔn)非冗余的基因數(shù)據(jù)庫(kù).(It provides the Sequences Record and Baseline for biological studies and it's a non-redundant set of reference standards derived)
- 用readseq對(duì)gb數(shù)據(jù)轉(zhuǎn)換：
efetch -db=nuccore -format=gb -id=AF086833 > AF086833.gb
###轉(zhuǎn)gb到cds序列fasta序列
readseq -p -format=FASTA AF086833.gb ## complete genome seq
readseq -p -format=FASTA -feat=CDS AF086833.gb
readseq -p -format=GFF AF086833.gb
readseq -p -format=GFF -field=CDS AF086833.gb
2. FASTA
就是第一行以'>'開始的序列的各種信息，第二行具體atgc
3. FASTQ
二代測(cè)序數(shù)據(jù)存儲(chǔ)格式鼠哥，第一行'@'熟菲，第二行''序列，第三行'+或者序列信息''朴恳，第四行測(cè)序的數(shù)據(jù)質(zhì)量抄罕。
Charpter 7: 獲取序列How to get DATA
1. 通過(guò)LINUX系統(tǒng)直接從NCBI上下載數(shù)據(jù)
### 以genebank格式序列下載
efetch -db=nuccore -format=gb -id=AF086833 > AF086833.gb
### 以fatst格式下載序列
efetch -db=nuccore -format=fasta -id=AF086833 > AF086833.fasta
### 取序列的某一段
efetch -db=nuccore -format=fasta -id=AF086833 -seq_start=1 -seq_stop=100
## 互補(bǔ)序列
-strand 1
-strand 2
- 通過(guò)對(duì)文章里的project accession 號(hào)來(lái)下載數(shù)據(jù)，以PRJNA257197為例子**
esearch -db nucleotide -query PRJNA257197
esearch -db nucleotide -query PRJNA257197 | efetch -format=fasta > genomes.fa
-db=protein
2. 下載SRA數(shù)據(jù)庫(kù)數(shù)據(jù)
-
SRA數(shù)據(jù)庫(kù)(short read Archive):數(shù)據(jù)存儲(chǔ)的格式
- Bioproject：PRJN開頭于颖，研究項(xiàng)目的介紹
- BioSample：SRS/SRMN開頭呆贿，生物樣品來(lái)源的介紹
- SRA Experiment: 單獨(dú)樣本的測(cè)序文庫(kù)
- SRR Run：直接相關(guān)數(shù)據(jù)存儲(chǔ)格式
==根據(jù)eseach 和 efetch直接下載數(shù)據(jù)== （都不用各種網(wǎng)站點(diǎn)了，實(shí)用下載數(shù)據(jù)新技能；衅Ｕケ馈！）
esearch -db sra -query SRP020911 | efetch -format runinfo >SRP020911.csv
## 查看csv第一行標(biāo)題信息
head -1 SRP020911.csv |tr ',' '\n'|nl ### 第一列為 SRR信息號(hào)章母，第十列為SRR的下載地址
### 獲取SRR的信息號(hào),SRR下載地址
cat SRP020911.csv| tail -n +2 | cut -d, -f 1
tail -n +2 | cut -d, -f 10
### 結(jié)合xargs直接完成數(shù)據(jù)的下載：
cat SRP020911.csv| tail -n +2 | cut -d, -f 1 | xargs -i prefetch {}
cat SRP020911.csv| tail -n +2 | cut -d, -f 10 | xargs -i echo wget {} \&
3. 對(duì)Seqkit工具的應(yīng)用
- 對(duì)fq文件基本情況統(tǒng)計(jì)
seqkit stat *.gz - 統(tǒng)計(jì)每個(gè)序列的GC情況
seqkit fx2tab --name --only-id --gc *.gz - ==根據(jù)id號(hào)提取序列==
### 隨機(jī)創(chuàng)造個(gè)id列表
seqkit sample -p 0.001 duplicated-reads.fq.gz |seqkit seq --name --only-id > id.txt
###根據(jù)id列表提取序列
seqkit grep --pattern-file id.txt duplicated-reads.fq.gz
或者-f參數(shù)
### 提取目標(biāo)序列
seqkit grep -p "AT1G30490.1" ath.pep.fa
- 去除掉含有簡(jiǎn)并堿基的序列如N/K等
### 生成含有簡(jiǎn)并堿基序列的id
seqkit fx2tab -n -i -a viral.1.1.genomic.fna.gz | csvtk -H -t grep -f 4 -r -i -p "[^ACGT]" | csvtk -H -t cut -f 1 > id.txt
### 去除id.txt母蛛，即反向匹配。
seqkit grep --pattern-file id.txt --invert-match viral.1.1.genomic.fna.gz > clean.fa
### 定位這些含有簡(jiǎn)并堿基的序列乳怎。
seqkit grep -f id2.txt viral.1.1.genomic.fna.gz |seqkit locate -i -P --pattern K+ --pattern N+
- 定位motif
## 以文件中的motif為標(biāo)準(zhǔn)定位
seqkit locate --degenerate --ignore-case --pattern-file enzymes.fa viral.1.1.genomic.fna.gz
### 根據(jù)自己的指定的motif序列定位彩郊，可包含簡(jiǎn)并堿基
seqkit locate -i -d -p CHWWWWWWDG seq.fa
- 根據(jù)序列的長(zhǎng)短從小到大排序
seqkit sort --by-length virus.fna.gz > virus.genomics.sorted.fa
文件過(guò)大用 --two-pass
Charpter 8：二代三代測(cè)序數(shù)據(jù)基本信息Sequencing instruments
DNA sequencing at 40: past, present and future
1. 測(cè)序儀器
- 一代測(cè)序：Sanger法，利用加不同熒光的雙脫氧核苷酸ddNTP會(huì)中止PCR蚪缀。
- 速度快秫逝，一次只能一條，1000bp~1500bp左右
- 二代測(cè)序：以Illumina為代表邊合成邊測(cè)序询枚。對(duì)隨機(jī)打斷成150~300bp的短片段 +上固定堿基的接頭(adapter)+標(biāo)簽(tag/index)违帆。
- 由于建庫(kù)中利用了PCR富集序列，因此有一些量少的序列無(wú)法被大量擴(kuò)增金蜀，造成一些信息的丟失刷后，且PCR中有概率會(huì)引入錯(cuò)配的堿基
- 三代測(cè)序：
- 單分子熒光測(cè)序(SMRT):
- 納米孔測(cè)序
- 無(wú)需擴(kuò)增，基于納米科技渊抄，對(duì)單分子鏈DNA/RNA直接用合成/降解/通過(guò)納米孔等方式直接測(cè)序
2. 錯(cuò)誤率
- ILLUMINA 0.1%的錯(cuò)誤率尝胆，理論上一條序列的測(cè)序準(zhǔn)確度是遞減的，也就是前面準(zhǔn)確度高护桦，后面低含衔。且在一些GC率高的區(qū)域，測(cè)序質(zhì)量會(huì)顯著降低二庵。
- PacBio 10%的錯(cuò)誤率
- MinION ：最高達(dá)到20%的錯(cuò)誤率贪染，且錯(cuò)誤多為固定錯(cuò)誤，無(wú)法通過(guò)測(cè)序深度來(lái)矯正催享。
可以看到三代測(cè)序的數(shù)據(jù)錯(cuò)誤率都很高抑进，所以常需要加測(cè)二代的數(shù)據(jù)來(lái)對(duì)三代數(shù)據(jù)矯正。
Charpter 9：測(cè)序質(zhì)控
分為兩步：第一步FastQC/MultiQC展示測(cè)序的質(zhì)量睡陪；第二步Trimmomatic/fastp去除低質(zhì)量序列
1. FastQC/MultiQC 顯示測(cè)序質(zhì)量
### 下載數(shù)據(jù)
wget http://data.biostarhandbook.com/data/sequencing-platform-data.tar.gz
tar xzvf sequencing-platform-data.tar.gz
### 直接fastqc顯示測(cè)序質(zhì)量
ls *.fq.gz |parallel fastqc -o ./ --nogroup {} &
### MultiQC批量顯示測(cè)序的質(zhì)量
##在已經(jīng)用fastqc得到測(cè)序的html文件和zip文件夾
multiqc *fastqc.zip --pdf
在運(yùn)行過(guò)程中曾出現(xiàn)too many levels 錯(cuò)誤寺渗，需要調(diào)用fastqc軟件的絕對(duì)地址來(lái)運(yùn)行程序
2. Trimmomatic/fastp質(zhì)量控制
- 用Trimmomatic 從5‘端開始切除低質(zhì)量堿基（LEADING：5）TRAILING:20表示切除堿基質(zhì)量小于20的堿基，MINLEN:50表示過(guò)濾切除后長(zhǎng)度小于50個(gè)堿基的reads
- 滑動(dòng)窗口：SLIDINGWINDOW:3：20表示設(shè)置窗口大小為3個(gè)堿基兰迫，從5'端開始切除reads信殊，read的平均質(zhì)量低于20，開始切除后邊的堿基序列汁果。過(guò)濾長(zhǎng)度小于50堿基的reads涡拘，使用更大的窗口可以減輕切除的力度
trimmomatic PE SRR11111_1.fq SRR11111_2.fq trimmed_1.fq unpaired_1.fq trimmed_2.fq unpaired_2.fq ILLUMINACLIP:/home/wangtianpeng/anaconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:5 MINLEN:20
trimmomatic PE SRR1553607_1.fastq SRR1553607_2.fastq trimmed_1.fq unpaired_1.fq trimmed_2.fq unpaired_2.fq SLIDINGWINDOW:4:30
- 去除特定的adapter
cat ~/miniconda3/envs/bioinfo/opt/fastqc-0.11.5/Configuration/adapter_list.txt
echo ">nextera" > nextera.fa
echo "CTGTCTCTTATACACATCTCCGAGCCCACGAGAC" >> nextera.fa
## 最后加上ILLUMINACLIP:參數(shù)
ILLUMINACLIP:adapter.fa:2:30:5
- 一鍵化軟件fastp
Github上地址https://github.com/OpenGene/fastp
###一鍵化輸出
fastp -i in.R1.fq -o out.R1.fq -I in.R2.fq -O out.R2.fq
### 默認(rèn)自動(dòng)化去接頭，如果不需要參數(shù)-A
### PE 數(shù)據(jù)的堿基校正 PE 數(shù)據(jù)的每一對(duì) read 進(jìn)行分析据德，查找它們的 overlap 區(qū)間鳄乏，然后對(duì)于 overlap 區(qū)間中不一致的堿基跷车，如果發(fā)現(xiàn)其中一個(gè)質(zhì)量非常高，而另一個(gè)非常低橱野，則可以將非常低質(zhì)量的堿基改為相應(yīng)的非常高質(zhì)量值的堿基值
-c 參數(shù)使用
fastp -q 30 -5 -l 100 -i il_1.fq.gz -I il_2.fq.gz -o i1_clean_1.fq -O i1_clean_2.fq
這里標(biāo)準(zhǔn)為：平均質(zhì)量高于Q30朽缴，對(duì)5‘端進(jìn)行低質(zhì)量堿基刪除，保留大于100bp的短讀水援。
fastp -3 -W 4 -M 20 -a AGATCGGAAGAG -i $i.fastq.gz -o ../clean_data/fastp_out_2.0/$i.qc.fq.gz
對(duì)3'端低質(zhì)量堿基刪除密强，-W -M 為sliding window參數(shù) -a 指定接頭序列
