要处理四个矩阵文件，分别为GSE109555 GSE136447 GSE36552 GSE66507
从ncbi网站上下载原始.sar文件，以及转化压缩成为.fastq.gz文件
在后台处理利用nohup&或者screen
cat GSE109555.txt | while read i;
do
/data/biosoft/software/sratoolkit.2.9.6-1-ubuntu64/bin/prefetch -X 9999999999999999999 -O `pwd` $i && echo "**${i}.sra done**";
fastq-dump --split-3 ${i}.sra;
gzip ${i}.fastq
done
————————————生成count矩阵，要求比对率为60%以上————————————
GSE66507，是双端测序文件，长度为50bp
比对过程中出现问题
多是由于数据没有下载完整造成的，此时需要重新下载
先是利用组里的index进行比对；结果还可以？大致有2/3比对率是大于百分之50的
GSE36552，是单端测序文件，长度为100bp
准确率较好，基本可以稳定在55%以上
但是对于
GSE109555，是双端测序文件，长度为150bp
以及
GSE136447，是双端测序文件，长度为150bp
比对率较差，都只有15%左右
更换比对文件Homo_sapiens.GRCh38.99.gtf;hg38_ercc.fa
STAR建立index文件
STAR --runMode genomeGenerate --genomeDir /data/shift/other/index_STAR --runThreadN 20 --genomeFastaFiles hg38_ercc.fa --sjdbGTFfile Homo_sapiens.GRCh38.99.2.gtf
出现问题：.fa和gtf文件里的chr1 和1之间有差异，求助解决
Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA fill
确认文件中的开头
grep -o -E "^\w+([-+.]\w+)*" Homo_sapiens.GRCh38.99.gtf|sort|uniq
将.gtf文件中加chr
继续建立index文件，而后进行比对：本次比对对数据进行了清洗
#!/bin/bash
cat /data/shift/other/GSE36552/GSE36552.txt|while read i;
do
fastp -i /data/shift/other/GSE36552/GSE36552_data/${i}.fastq.gz -o /data/shift/other/GSE36552/GSE36552_data_clean/${i}_clean.fastq.gz;
STAR --runThreadN 12 --genomeDir /data/shift/other/index_STAR --readFilesCommand zcat --readFilesIn /data/shift/other/GSE36552/GSE36552_data_clean/${i}_clean.fastq.gz --sjdbGTFfile /data/shift/other/index_STAR/Homo_sapiens.GRCh38.99.2.gtf --outFileNamePrefix /data/shift/other/GSE36552/GSE36552_count.38.99/${i}_ --outSAMtype BAM SortedByCoordinate;
featureCounts -p -t exon -g gene_id -a /data/shift/other/index_STAR/Homo_sapiens.GRCh38.99.2.gtf -o /data/shift/other/GSE36552/GSE36552_count.38.99/${i}_counts.txt /data/shift/other/GSE36552/GSE36552_count.38.99/${i}_Aligned.sortedByCoord.out.bam;
done
但是比对结果并不尽如人意，比对成功率略低于组里文件：思考建立STAR index时是否应该根据不同的序列长度增加参数 --sjdbOverhang reads长度的最大值减1，默认是100
提出解决办法：随机选取一对文件fastq从中取出50bp(两个文件对称选取)的长度，与人进行比对
先对GSE136447尝试
利用软件seqkit
从Download - SeqKit - Ultrafast FASTA/Q kit (shenwei.me)下载软件
操作方法指路：seqkit：序列梳理神器-统计、格式转换、长度筛选、质量值转换、翻译、反向互补、抽样、去重、滑窗、拆分等30项全能..._刘永鑫Adam的博客-CSDN博客
报错：seqkit无法处理.gz文件
报错：更换剪切长度50，100都报错
截取后文件分析
比对源文件
原文件就有点问题（心酸）
尝试解决途径，仅保留后50位数据，报错不行
更改保留前100数据尝试,比对成功
突然发现小时错了/(ㄒoㄒ)/~~看比对率的文件应该是，*_Log.final.out
GSE66507的比对率为80%以上
GSE36552的比对率为87%以上
现有问题GSE136447 featureCounts的assignment低
按照文章方式进行比对
利用hisat2
我知道为啥我的GSE109555一直都不好使了
原因一：有好多数据都没下完整（可能是由于网络原因），这部分需要重新下载
原因二：文章中说不能直接用star，需要将每一个文件拆分后，将子文件做star
1.利用 zcat SRR8569907_2.fastq.gz | wc -l 与原文件比对（一个一个比啊，愿天堂没有数据不完整）
2.下载 GSE109555_TrioSeq_DataInfo.txt.gz 作为输入执行
perl s01.Barcode_UMI_QC_per1w_V2.pl test_1.fq.gz test_2.fq.gz 不同的文件对应不同应当根据SRR对应的GME查找 ./
打开文件夹以后里面对应的是每个压缩文件
而后就可以执行star了！
85%还挺好，可以说是感动中国了/(ㄒoㄒ)/~~接下来就是等待~
网友评论